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ABSTRACT 

Using  results  of  DNS  in  the  case  of  two-dimensional  homogeneous  isotropic 
flows,  we  first  analyze  in  detail  the  behavior  of  the  small  and  large  scales  of 
Kolmogorov  like  flows  at  moderate  Reynolds  numbers.  We  derive  several 
estimates  on  the  time  variations  of  the  small  eddies  and  the  nonlinear  in¬ 
teraction  terms;  those  terms  play  the  role  of  the  Reynolds  stress  tensor  in 
the  case  of  LES.  Since  the  time  step  of  a  numerical  scheme  is  determined  as 
a  function  of  the  energy-containing  eddies  of  the  flow,  the  variations  of  the 
small  scales  and  of  the  nonlinear  inter2w;tion  terms  over  one  iteration  can  be¬ 
come  negligible  by  comparison  with  the  accuracy  of  the  computation.  Based 
on  this  remark,  we  propose  a  multilevel  scheme  which  treats  dilferently  the 
small  and  the  large  eddies.  Using  mathematical  developments,  we  derive  esti¬ 
mates  of  all  the  parameters  involved  in  the  algorithm,  which  then  becomes  a 
completely  self-adaptive  procedure.  Finally,  we  perform  realistic  simulations 
of  (Kolmorov  like)  flows  over  several  eddy-turnover  times.  The  results  are 
analyzed  in  detail  and  a  parametric  study  of  the  nonlinear  Galerkin  method 
is  performed. 


'  Part  of  the  work  was  done  while  the  second  and  third  authors  were  visiting  the  Institute 
for  Computer  Applications  in  Science  and  Engineering  (ICASE),  NASA  Langley  Research 
Center,  Hampton,  VA  23681. 


in 


Introduction 


A  turbulent  flow  can  be  characterized  by  a  spatial  and  temporal  chaotic  be¬ 
haviour.  Indeed,  strong  gradients  may  appear  and,  in  the  physical  space, 
several  structures  with  various  sizes  are  generated  by  the  external  force  and 
by  the  flow  itself.  The  structures  are  convected  by  the  mean  flow.  In  the  case 
of  viscous  incompressible  flows  in  dimension  two,  thin  and  elongated  sheets  of 
vorticity  appear.  These  structures  are  characteristic  of  incompressible  flows. 
The  phenomenological  theory  of  turbulence,  introduced  by  Kolmogorov  in 
dimension  three  and  Kraichnan  in  dimension  two,  predicts  that  the  size  of 
the  smallest  scales  of  a  flow  decreases  while  the  nondimensional  Reynolds 
number  grows.  So,  the  number  of  degrees  of  freedom  required  to  describe 
the  motion  can  be  estimated  in  terms  of  the  Reynolds  number  (~  cRe  in 
space  dimension  two). 

Hence,  a  Direct  Numerical  Simulation,  i.e.  the  resolution  of  all  physically 
relevant  scales,  can  not  be  achieved  when  the  Reynolds  number  becomes  too 
high  and  then,  when  the  turbulence  is  fully  developed.  A  model  is  then  used 
in  many  simulations  ;  i.e.  the  small  scales  are  not  exactly  integrated  but 
their  interaction  with  the  large  ones  are  taken  into  account  in  a  simplifled 
way.  Basically,  the  aim  of  these  models  is  to  recover  the  large  eddies  of  the 
flow  (or  their  statistic)  without  explicitly  computing  all  the  scales  of  motion. 

In  relation  with  recent  developments  in  the  theory  of  dynamical  systems 
and  its  application  to  turbulence  phenomena,  new  objects  have  been  in¬ 
troduced  (exact  or  approximate  inertial  manifolds,  see  Foias,  Manley  and 
Temam  [1])  and  new  numerical  methods  have  been  proposed  (the  nonlinear 
Galerkin  methods,  see  Marion  and  Temam  [2]  and  [3],  the  incremental  un¬ 
known  method,  see  Temam  [4]).  These  objects  and  methods  are  based  on  a 
decomposition  of  the  unknown,  such  as  the  velocity  field,  into  its  small  scale 
component  and  its  large  scale  one,  and  essential  in  the  nonlinear  Galerkin 
method  is  a  systematically  differentiated  treatment  of  the  small  and  large 
scales. 

These  approximate  inertial  manifolds  are  subsets  of  the  phase  space  and 
consist  in  an  approximate  form  of  the  small  scale  equations.  They  provide 
an  adiabatic  law  modeling  the  interaction  between  the  low  and  the  high  fre- 
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quency  components  of  the  flow  ;  the  small  eddies  are  in  fact  expressed  as 
a  nonlinear  function  of  the  large  ones,  they  are  slaved  by  the  large  ones. 
Moreover,  these  manifolds  enjoy  the  property  that  they  attract  all  the  orbits 
exponentially  fast  in  time  and  that  they  contain  the  attractor  in  a  thin  neigh¬ 
borhood.  In  that  sense,  they  provide  a  good  way  to  approximate  the  solutions 
of  the  Navier-Stokes  equations.  The  nonlinear  Galerkin  method,  proposed 
by  Marion  and  Temam  [2],  consists  in  looking  for  a  solution  lying  on  these 
specific  subsets  of  the  phase  space.  Several  implementations  of  this  scheme 
have  previously  been  done  by  Jauberteau,  Rosier  and  Temam  ([5]  and  [6]), 
and  Dubois,  Jauberteau  and  Temam  ([7]).  Results  presented  in  these  papers 
were  mainly  feasibility  tests  involving  exact  solutions  (analytical),  where  the 
authors  tried  to  recover  known  velocity  and  pressure  fields.  These  tests  pro¬ 
vided  good  indications  on  computational  feasibility  and  performances,  but 
had  no  physical  relevance.  In  the  present  article  we  intend  to  develop  further 
the  study  of  the  numerical  implementation  of  the  nonlinear  Galerkin  method 
for  flow  problems.  We  have  several  objectives  which  we  describe  hereafter 
in  some  details.  Firstly,  we  would  like  to  compute  physically  more  relevant 
flows,  namely  here  Kolmogorov  type  flows.  Secondly,  the  effective  implemen¬ 
tation  of  the  nonlinear  Galerkin  method  involves  several  parameters  (such 
as  the  cut-off^  wavelength  between  low  and  high  frequencies,  and  the  time 
during  which  the  high  frequencies  are  allowed  to  be  frozen)  ;  and  another 
of  our  aims  is  to  conduct  a  parametric  study  of  the  method,  and  to  develop 
simple  self-adaptive  methods  for  the  determination  of  these  parameters. 

Here,  we  consider  two-dimensional  periodic  flows  governed  by  the  incom¬ 
pressible  Navier-Stokes  equations.  Of  course,  such  flows  correspond  more  to 
an  idealized  situation  rather  than  a  physical  one.  Nevertheless,  this  model 
contains  several  difficulties  which  occur  when  simulating  turbulent  flows,  and 
then  it  provides  a  good  test  to  implement  a  numerical  algorithm.  More  phys¬ 
ically  relevant  situations,  such  as  three-dimensional  flows,  will  be  presented 
elsewhere.  In  a  parallel  effort,  we  are  also  treating  the  case  of  a  flow  driven 
by  a  constant  pressure  gradient  between  two  infinite  parallel  plates,  namely 
the  channel  flow  (  see  for  instance  [8]). 

In  the  first  section  of  the  paper,  we  introduce  the  equations  and  several  nota¬ 
tions.  Then,  we  define  the  decomposition  of  the  velocity  field  into  large  scale 
components  and  small  scale  structures.  This  decomposition  depends  on  a 
cut-off  value,  corresponding  to  a  wave-number  of  the  Fourier  decomposition 
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of  the  unknowns,  and  provides  two  grids  in  the  physical  space.  In  previous 
theoretical  works  (see  for  instance  [1]  and  [9]),  the  authors  have  rigorously 
proved  that,  for  a  sufficiently  large  wavenumber,  the  large  scale  components 
contain  most  of  the  energy  of  the  flow.  The  small  scale  structures  repre¬ 
sent  then  small  quantities,  which  have  to  be  taken  into  account  in  order  to 
correctly  compute  the  large  scales  of  motion.  After  recalling  these  results 
and  comparing  them  with  numerical  simulations,  we  show,  based  on  both 
theoretical  and  numerical  inve.stigations,  that  the  time  variation  of  the  small 
eddies  over  one  time  step  can  be  smaller  than  the  energy  carried  by  these 
scales  themselves,  when  the  cut-off  value  is  chosen  sufficiently  large.  A  sim¬ 
ilar  result  is  also  proved  for  the  nonlinear  interaction  terms  expressing  the 
action  of  the  small  eddies  over  the  large  ones.  Then,  introducing  the  expected 
accuracy  of  the  computation  £  as  a  parameter,  we  deduce  that  there  exists 
a  level  such  that  the  one-step  time  variation  of  the  small  scale  components 
is  smaller  than  e.  We  note  here  that  this  result  is  not  in  contradiction  with 
the  fact  that  the  small  eddies  evolve  faster  than  the  large  ones.  In  fact,  their 
time  scales  are  smaller  but  their  order  of  magnitude  also  decreases  for  large 
wavenumbers. 

Hence,  the  time  variations  of  part  of  the  spectrum  can  be  locally  neglected. 
We  have  also  proved  that  the  interaction  terms  enjoy  the  same  property. 
These  two  fundamental  results  are  one  of  the  keys  of  the  numerical  algo¬ 
rithm  proposed  in  the  following  sections. 

Based  on  these  results,  we  have  implemented  a  spatial  and  temporal  mul¬ 
tilevel  adaptative  method  ;  i.e.  the  level  of  refinement,  which  define  the 
separation  of  the  flow  into  low  and  high  frequencies,  evolves  in  time  under¬ 
going  successive  multilevel  V-cycles  ;  moreover,  it  follows  the  dynamic  of  the 
small  and  the  large  scales  of  the  motion.  Moreover,  some  part  of  the  spec¬ 
trum  are  frozen,  as  well  as  the  interaction  terms,  during  short  time  periods. 
After  a  complete  description  of  this  multilevel  procedure,  we  derive  several 
estimates  of  characteristic  time  scales  for  the  small  structures  and  for  the 
terms  involving  their  interactions  with  the  large  scale  components.  These 
estimates  provide  an  efficient  way  to  evaluate  the  length  of  the  time  interval 
during  which  the  small  scales,  corresponding  to  several  levels  of  refinement, 
as  well  as  the  corresponding  interaction  terms  are  allowed  to  be  frozen  with¬ 
out  introducing  an  error  larger  than  the  accuracy  e.  Hence,  we  obtain  a 
completely  self-adaptative  procedure  which  depend  on  only  one  parameter. 
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namely  e.  At  the  end  of  each  frozen  period,  the  velocity  field  is  projected 
on  an  Approximate  Inertial  Manifold  (AIM),  which  provides  a  simple  find 
efficient  way  to  directly  compute  the  high  frequencies.  We  analyze  here  the 
perturbation  introduced  by  using  first  order  Approximate  Inertial  Manifolds 
and  we  deduce  in.  which  range  of  the  spectrum  the  manifold  can  be  used  to 
compute  the  small  scales. 

The  last  section  is  devoted  to  the  description  and  the  analysis  of  several 
numerical  simulations  performed  with  the  multiscale  (Nonlinear  Galerkin) 
method.  The  integral  scale  Reynolds  number  Rti,  vias  successively  set  to 
784  and  to  6,328  for  spatial  resolutions  of  respectively  (256)*  and  (512)*. 
In  these  cases,  a  full  dissipation  range  is  resolved  as  the  larger  wavenumber 
is  at  least  6  times  larger  than  the  dissipation  wavenumber.  The  Reynolds 
numbers  are  not  sufficiently  large  for  the  solutions  to  have  a  k~^  energy 
spectrum  power  range  ;  indeed,  the  computed  velocities  have  an  intermedi¬ 
ate  energy  spectrum  between  k~*  and  k~^.  Nevertheless,  the  different  scale 
components  are  strongly  time  dependent,  and  they  provide  interesting  tests 
for  the  multilevel  adaptative  procedure.  More  physically  relevant  situations, 
eg  simulations  with  larger  Reynolds  number,  will  be  treated  and  presented 
elsewhere.  For  each  simulation,  the  code  has  run  during  more  than  50, 000 
time  iterations,  which  represents  10s  of  unit  hours  on  a  Cray-2. 

In  this  last  section,  we  first  compare  the  results  obtained  with  a  previous 
version  of  the  algorithm  with  those  obtained  with  the  multiscale  procedure 
described  in  subsection  2.4.  For  this  comparison  a  Direct  Numerical  Sim¬ 
ulation  performed  with  a  standard  pseudo-spectral  Galerkin  approximation 
is  used  as  a  reference.  In  the  previous  versions  of  the  algorithm  (  [7],  [10]), 
there  was  no  control  of  the  variations  of  the  small  eddies.  Hence,  the  cut-off 
wavelengths  as  well  as  the  length  of  the  period  during  which  the  variations 
of  the  small  scales  are  frozen  were  not  properly  evaluated.  This  led  to  an 
accumulation  of  perturbations  and  a  loss  of  the  accuracy  ;  this  is  no  more 
the  case. 

Finally,  in  order  to  study  the  effect  of  the  cut-off  wavelength  on  the  numeri¬ 
cal  results,  we  perform  several  simulations  with  decreasing  cut-off  values.  In 
fact,  we  show  that  the  separation  wavenumber  can  be  taken  of  the  order  of 
the  dissipation  wavenumber  A:,,  ;  as  A:,,  is  of  the  order  of  20,  in  our  computa¬ 
tion,  the  cut-off  value  can  not  reasonably  take  smaller  value.  Nevertheless, 
the  standard  Galerkin  pseudo-spectral  method,  with  a  larger  wavelength  of 
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the  order  of  the  dissipation  one,  does  not  provide  a  sufficient  resolution  ;  the 
large  scales  of  motion  are  not  recovered  in  that  case.  Hence,  the  multiscale 
method  provides  an  efficient  way  to  simulate  the  enstrophy  dissipation  and 
the  energy  dissipation  ranges.  In  order  to  describe  the  effect  of  the  multiscale 
procedure  inside  the  power  range,  simulations  with  larger  Reynolds  number 
have  to  be  performed.  Such  study  will  be  presented  elsewhere. 
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1  Motivations. 

1.1  Preliminary. 

We  consider  here  two-dimensional  viscous  incompressible  flows  driven  by  an 
external  volume  force  f ,  and  governed  by  the  Navier-Stokes  equations  : 

’  ^-j/Au  +  (u- V)u  + Vp  =  f, 

•  V-u  =  0,  (1) 

,  u(x,t  =  0)  =  Uo(x), 

where  u(x,<)  =  (ui(x,<),U2(x,<))  is  the  velocity  field,  p(x,<)  the  pressure 
(x  =  (xi,a;2)),  u  is  the  kinetic  viscosity  and  the  density  is  set  to  unity. 
Ek]uation  (1)  is  supplemented  with  boundary  conditions,  namely  u  and  p  are 
periodic  of  period  Li  (resp.  L2)  in  the  direction  Xj  (resp.  x^).  We  denote 
by  fi  =  (0,  Li)  X  (0,  L2)  the  period  and  we  assume  that  the  average  of  the 
external  force  over  the  period  ft  vanishes,  i.e.  : 

=  0,  Vt. 

Taking  the  average  over  the  whole  domain  ft  of  the  momentum  equation 
in  (1)  and  using  the  periodicity  of  u  and  p,  we  obtain  : 

Assuming  now  that  the  average  of  the  initial  condition  is  also  zero,  we  con¬ 
clude  that  the  average  of  the  velocity  field  vanishes  : 

=  0,  V(. 

We  now  denote  by  u;  =  V  x  u  the  vorticity  and  as  usual  we  rewrite  the 
conservation  of  momentum  equation  in  (1)  as  : 

-^-i/Au-b(u;xu)-l-VP  =  f,  (2) 
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where  P  =  p  +  ^  I  u  p  .  This  choice  is  motivated  by  two  facts.  First,  it 
is  important  to  note  that  a  pseudo-spectral  evaluation  of  (a;  x  u)  requires 
three  Fast  Fourier  Transforms  (FFT)  fewer  than  the  usual  form  (u  •  V)u. 
Also,  this  form  of  the  equations  semi-couserves  the  kinetic  energy  for  iuviscid 
flows,  when  a  collocation  or  pseudo-spectral  discretization  is  used  (see  [11]) 
and  ensures  numerical  stability,  while  the  original  form  does  not. 


Since  the  unknowns  u  and  p  are  space  periodic  functions,  they  can  be 
expanded  in  Fourier  series  : 


u(x,t)  =  ^  u(k,0  tnk(x), 
k€Z^ 

p(x,<)  =  p{k,  t)  ti;k(x), 
keZ* 


(■■i) 


where  iwic(x)  =  =  {k\l Lx,k2l L2)  and  ki,  x  =  {k\Xx/Li  +  k2X2l L2)  • 

The  Fourier  coefficients  u(k,  t)  =  (resp.  p(k,  t))  are  com¬ 

plex  numbers  and  we  recall  that  : 


Ui(-k,t)  =  ?  =  1  or  2, 

p(-k,0  =  p(k,t), 


(4) 


where  c  denotes  the  complex  conjugate  of  c.  This  property  is  very  important 
in  practice  ;  it  allows  to  store  twice  fewer  coefficients  in  the  second  direction 
of  the  spectral  space,  when  we  look  for  a  finite  approximation  of  u(x,  t)  and 
p(x,  t),  i.e.  by  storing  u(k,  t)  for  k  such  that  k2  >  0,  (4)  allows  us  to  recover 
the  Fourier  coefficients  u(k,  t)  for  kj  <  0. 


Let  us  now  introduce  Fdiv  the  projection  operator  onto  the  divergence 
free  space  ;  P^y  can  be  easily  expressed  as  : 

Pdiy'Pix)  =  5^  (<^k  -  7^(k  •  <^k))  «nc(x). 

k€Z^  V  I  k  I  / 

Assuming  that  u  and  p  lie  in  the  proper  Hilbert  spaces  (see  for  instance  [12],  [13]) 
and  applying  Penv  to  the  Navier-Stokes  equations,  we  obtain 

^-uAu  +  B{u,u)  =  g,  (5) 
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where  g  =  Pdivf  B{u,  u)  is  a  bilinear  form  defined  by  : 

B(u,  u)  =  Pdiv(u»  X  u) 

=  S  ((w'^u)k-7^(k  (w'^u)k))  u;k(x)  ^  ^ 
kgZ*  I  1  >' 

As  we  shall  see  later,  the  numerical  procedures  are  directly  applied  to  this 
last  form  of  the  Navier-Stokes  equations. 

1.2  Separation  of  scales. 

1.2.1  Small  and  large  scales  equations. 

Let  us  choose  an  integer  N.  We  introduce  the  orthogonal  projector  onto 
the  space  spanned  by  the  first  Fourier  modes.  A  pseudo-spectral  Galerkin 
approximation  ua^  of  the  velocity  field  u  is  given  by 

u/v(x,t)  =  Y,  u(k,0»%(x), 
k€/w 

where  /yv  =  [1  —  N/2,  N/2]  x  [0,  N/2] ,  and  satisfies  the  following  equation 

-  i/Auyv  +  PyvP(uyv>ujv)  =  Pyvg-  (7) 

Here  we  only  consider  external  force  with  no  high  wavenumber  coefficients, 
namely  : 

g(k)  =  0,  for  k  €  /Ar\/mo»  for  some  mo.  (8) 

Instead  of  (8),  we  can  also  assume  that  the  spectrum  of  the  externa]  force 
has  an  exponential  decay  for  wavenumbers  larger  than  a  given  one. 

We  now  introduce  the  decomposition 

uat  =  yyv,  +  ZN^,  with  mo  <  Ni  <  N,  (9) 

where 

YNi  =  Pat,  Uat  =  u(k,t)ti;k  (10) 

ke/w 
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and 

u(k,0«^-  (11) 

We  note  that,  in  the  decomposition  (9),  ysi  represents  the  large-scale  struc¬ 
tures  of  the  flow,  and  Zjv,  corresponds  to  the  small-scale  components,  or  in 
other  words,  the  fine  structures  of  the  flow.  Since,  the  operators  Pat,  and 
Qjvj  commute  with  the  differentiation  operators,  we  have  : 

P/v,(Au^)  =  A(Pv,UAf)  =  Ayyvi, 

and 

(Au;v)  =  A(Qj^jUjv)  =  AZa/j  . 

By  projecting  (7)  with  respect  to  Psi  and  we  obtain  the  following 
coupled  system  of  ordinary  differential  equations  (ODE)  : 

-  vAyst  +  B{yNi  +  ^Ni^YNi  + 

(12) 

^  ZNj,yw, +ZN,)  =  0. 

The  initial  conditions  associated  with  (12)  are  Pv,Uo  and  Qjv,Uo.  In  (12),  as 
well  as  in  the  following,  we  will  only  use  the  projection  operator  Pnj  with 
N-i  larger  than  mo- 

We  now  introduce  two  norms  : 

\<fi\2=  ^1  y>(x,0  r  dx)  , 

and 

(/jv».(x,01’  *t)'', 

for  any  given  field  ^(x,  t)  —  (v?i(x,  <),  ¥>2(x,  <)).  We  note  that  these  norms  are 
related  to  well-known  physical  quantities.  Indeed,  with  the  above  definitions, 
the  kinetic  energy  e(u)  of  the  velocity  field  u  is  related  to  |  —  I2  by 

e(u)  = 
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and  in  the  case  considered  here  of  periodic  boundary  conditions,  the  enstro- 
phy,  ens(u)  is  : 

en5(u)  =  ^11  u  11^  =  ^  1  u>  |:j  . 

We  now  recall  theoretical  results,  constituting  the  basis  of  this  framework  ; 
they  were  established  by  Foias,  Manley  and  Temam  in  [1]  and  [9].  In  [1],  the 
authors  proved  that  for  sufficiently  large  values  of  N  and  Ni ,  namely  N  and 

Ni  ~  where  G  designates  the  Grashoff  number  G  =  (Ai  ~  \/L^ 

where  L  is  a  characteristic  length  of  the  domain),  and  after  a  transient  time, 
depending  only  on  the  data,  the  small  scale  components  (t)  remain  small 
in  both  norms  introduced  above,  while  the  large  scales  yN^  {t)  are  of  the  order 
of  u(<).  We  do  not  know  how  close  to  the  inertial  range  the  corresponding 
cut-off  eigenvalue  Ayv,  =  cN^  can  be,  and  it  is  one  of  our  objectives  in  this 
work  to  explore  this  point.  In  [9],  the  authors  derived  other  estimates  of  the 
ratios 

I  ZNi(0  I2  ,  IIZAf.W  II 
|yAr,(0l2  Ilyvv.WII 


in  terms  of  the  Grashoff  number  at  appropriate  powers  of  1/ (Ni +1 ),  when  the 
cut-off  value  Ni  is  still  of  the  order  of  G^.  Moreover,  by  using  the  analyticity 
in  time  of  the  solutions  of  the  Navier-Stokes  equations  and  invoking  the 

Cauchy  formula,  one  can  show  that  |  |2>  as  well  as  ||  ||>  are 

respectively  of  the  order  of  |  z^r,  I2  and  ||  Za?,  ||,  for  Ni  sufficiently  large, 
corresponding  to  a  wavenumber  in  the  dissipation  range.  Even  if  these  results 
do  not  provide  fine  estimates  on  the  norms  of  the  time  derivative  of  the  small- 
scale  components,  they  show  that  the  time  derivatives  of  the  velocity  field  u 
have  also  a  decaying  spectrum,  which  is  not  the  case  of  course  in  the  inertial 


range. 

This  behavior  of  the  small-scale  z^r,  is  very  important  and  is  one  of  the 
keys  of  the  numerical  implementation  of  the  multilevel  method,  Jis  we  will 
see  in  the  following  sections.  Although,  these  results  were  proved  theoreti¬ 
cally  only  when  Ni  is  of  the  order  of  G^,  and  Ajv,  is  far  in  the  dissipation 
range,  numerical  results  show  that  za^,  is  small  compared  to  yA/,  for  much 
smaller  values  of  Ni,  and  the  same  is  true  for  their  time  derivative.  From 
the  phenomenological  theory  of  two-dimensional  turbulence  point  of  view, 
due  in  part  to  Kraichman  (cf  [14]  and  [15]),  the  energy  spectrum  of  a  de- 


10 


veloped  turbulent  flow  is  expected  to  display  a  power-law  inertial  range  and 
an  exponential  decay  for  wavenumbers  larger  than  the  Kraichnau’s  dissipa¬ 
tion  wavenumber  (see  Figure  1).  So,  the  theoretical  results  mentioned  above 
are  in  agreement  with  this  assumption  in  the  sense  that  for  sufSciently  large 
wavenumbers,  the  energy  spectrum  has  a  decaying  behaviour.  The  exponen¬ 
tial  range  is  generally  referred  as  the  dissipation  range. 


Figure  1:  Energy  Spectrum  E{k). 

Due  to  the  bilinearity  of  B,  we  can  split  B{un,Un)  into 
B(uiv,  uyv)  =  B{yisf,  +  zjv, ,  y/v,  +  ) 

=  B(y«, .  yw, )  +  BtaiCyw,  ,*«,), 


(13) 


where 


BintiVNi  ,^Ni)  =  5(y^, ,  ZiV, )  +  B{zNi  ,  y^, )  +  ,  Z;v,  )• 

It  represents  the  interaction  between  the  small  and  the  large  structures  and 
the  interaction  of  the  small  structures  among  them.  We  now  recall  that  a 


11 


classical  (Galerkin)  approximation  of  u  based  on  the  first  Nf  Fourier  coef¬ 
ficients  consists  in  setting  Zni  to  zero  in  (12),  so  that  PNiBaaiySit^Ni)  is 
neglected.  If  zsi  represents  a  small  part  of  the  kinetic  energy,  namely 

c(z/v,)  <  c(yAf,), 

then,  we  can  expect  that 

I  PSiBintiySi^^Ni)  U  <■  I  PNiB{yNi,yNi)  b  • 

Therefore,  the  contribution  of  the  interaction  terms  is  much  smaller  than  the 
component  of  the  nonlinear  terms,  involving  only  the  large  scales  This 
argument  can  be  justified  when  comparing  Figures  4  and  5.  More  precisely, 
we  have  the  following 

I  Pst  Bmt(y7v, ,  a/v, )  b  <  I  Pni  B(yjv, ,  y^, )  b 


I  Ay/v,  b 

However,  the  interactions  terms  can  not  be  neglected  in  the  first  equation  (12) 
imless  Ni  is  taken  very  large.  Indeed,  their  effects  on  long  time  computations 
can  modify  the  behavior  of  the  large  scale  components.  Since  the  evaluation 
of  these  terms  at  each  time  step  is  very  expansive,  we  try  to  take  their  effects 
into  account  in  a  simplified  manner.  We  will  see  that  the  time  variations 
of  zni  and  Pni  5int(yAf, ,  Zyv, )  are  locally  negligible,  and  this  will  allow  us  to 
freeze  them  during  some  parts  of  the  iterations. 


1.2.2  Time  variations  of  z/v,  and  Pi<fiBiat{yNi,ZNi)- 

We  now  aim  to  derive  an  estimate  of  the  variations  of  Za^,  over  one  time 
iteration.  This  quantity  can  be  represented  by 

Ayv,  =  At  I  ZN^  b, 

where  the  dot  represents  the  differentiation  with  respect  to  t.  The  time  step 
At  is  given  by  the  CFL  stability  condition  : 

At  N  \  un  |l«><  with  a  <  1.  (14) 
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We  recall  that  N  is  the  total  number  of  modes  in  each  direction.  We  assume 
that  the  smallest  scale  Is  =  ^jks-,  ks  =  N/2,  is  smaller  than  the  Kraichnan 
dissipation  scale  t^,  =  I /k„.  Let  us  first  consider  the  case  where  ts^  =  1  /ksi  is 
lying  between  if,  and  tsj  f-Ni »  f-n  are  of  the  same  order  of  magnitude  : 

f-S  <  ^Ni  £ 

Let  us  assume,  as  shown  in  Figure  6,  that  the  time  derivative  zsi  is  of  the 
order  of  the  dissipation  term  (remember  we  are  in  the  dissipation  range), 


I  I2  ~  1/ 1  Azjv,  I2 . 

Due  to  the  exponential  decay  of  the  velocity  spectrum  in  the  dissipation 
range,  we  can  write 

V  I  Azn,  I2  <  Cl  vk\  I  zs,  I2,  (15) 

where  ci  is  a  nondimensional  constant  of  the  order  of  unity.  It  follows  by  (14) 
that  : 

=  At  I  ZNi  I2  <  C2  I/A<  I  Zjv,  I2 

(16) 

Then,  using  the  estimate  of  k„  obtained  in  [16]  namely 

-  ^2/3aJ/3  ’ 


we  obtain  : 


A  ^  I  8  l2^  (  Hli  \  1/3  I  _  I 

An,  <  C2  , - ; —  I  -r-*-  ]v  '  |  Zn,  |a 

\/2Ai  I  lioo  \ksknJ 


As  ksi  <  ks,  we  obtain  ; 


A  ^  1  8  l2^  1-*^®  I  T 


Since  fcjv,  ~  fc,,  then  for  sufficiently  small  values  of  the  viscosity  v,  the  vari¬ 
ations  of  zjv,  over  one  time  iteration  are  much  smaller  than  |  Zjv,  I2  • 
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We  now  try  to  derive  a  similar  estimate  when 


—  ^Si  > 

i.e.  when  tsi  is  in  the  inertial  range.  In  that  case,  we  assume  as  shown 
on  Figure  7,  that  the  time  derivative  i/Vj  is  of  the  order  of  the  interaction 
nonlinear  terms 

I  iSi  b  ~  1  ^jv,^int(yN,,ZjV,)  I2  . 

We  recall  that 

I  Qs^Bintiysn^Ni)  I2 

=  I  Qsi B{yNt ,^Ni)  +  Qn, B{zn, , yjv, )  +  Qn, B{zn, , zjv, )  I2 

(17) 

^  I  QNiB{yNi,ZNi)  I2  +  I  QNiB{ZNi,ySr)  I2 
+  I  QNtB{^Ni,ZNi)  I2  • 

As  is  a  projection  operator,  we  have 

I  QN^B{yN^,ZN^)  I2  <  I  B{yNi,Zs,)  is  . 

The  bilinear  form  B  can  be  estimated  as  follows  (see  for  instance  [12]  and  [13]) 

I  Biys^zs,)  I2  <  C4  1  ys,  |l~||  zn,  || . 

As  I  yNi  Uiv  \l  00  j  for  sufficiently  large  values  of  N\ ,  we  obtain 

I  B{yN^,ZNi)  I2  <  C4  I  |i,oo||  zat,  || . 

The  decay  of  the  velocity  Fourier  components  implies  that 
IjaN.li  <  Ciksi  I  ZN,  12,  where  C5  ~  1, 

so  that 

I  ^jv,^(y/V„Z/V,)  I2  <  C4  C5  ksi  I  ujv  |l<»|  Ziv,  (2  . 

We  also  have 

I  QNiB{ZNx,yNx)  I2  <  I  ^(ZAf,,y7V,)  I2, 
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which  can  be  estimated  by  (see  [12],  [13]) 

I  B{zNi,yNi)  b  <  C6  I  zjv,  jj''*]  Az^,  \y'^  l|yjv,|| . 

In  fact,  the  norm  jj  ysi  ||  and  the  infinity  norm  |  y^v,  l^oo  are  of  the  same 
order.  We  then  obtain 

I  B{zNi,yNi)  b  <  C7  I  zsi  Azjv,  ysi  |l«  ■ 

From  the  inequality  (15),  we  deduce  that 

1  QN,^(*JVi,yNj)  b  ^  cg  fcw,  I  ysx  [l"!  zjvi  I2  • 

A  similar  estimate  can  be  derived  for  the  third  term  Q’^^B{zsxi^Nx)-  We 
finally  obtsun 

I  Qjv,  Biotiysi ,  zjv, )  b  <  C9  k^i  I  utv  |l~  |  zn,  b  * 


Then, 

Ayv,  <  Cg  AryVjAt  I  Uyv  |£,~|  ZjV,  b  • 

As 

^  a  _  Q _ 

~  I  Uyv  |l~  N  I  Un  |l“>  \/2  ks  ’ 

we  find 

^  {^) ' 

which  also  implies  that  Ayvi  <|  ZyV|  b  inertial  range  as  well.  Fig¬ 

ure  8  supports  the  theoretical  estimates  derived  above.  We  want  to  note, 
at  this  point,  that  these  results  will  remain  valid  if  energy  is  injected  in  the 
small  scales  by  a  given  external  force,  having  a  decaying  spectrum  ;  i.e.  if 
g(k)  I  k  )“  for  I  k  I  larger  than  a  given  wavenumber. 


Let  us  now  introduce  e  the  accuracy  of  the  computations  ;  e  represents 
here  an  energy  error.  We  assume  that  e  is  a  given  parameter,  which  can 
be  chosen  under  several  considerations  ;  examples  will  be  given  in  Section  3 
devoted  to  the  description  of  the  numerical  results.  The  accuracy  e  can  be 
associated  with  a  small  scale,  i.e. 

there  exists  a  level  Nx{e)  (<  N)  such  that  . 
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According  to  the  estimate  previously  derived,  ^Ni{e)  is  much  smaller  than  e. 
From  the  previous  estimates  we  can  deduce  that  there  exists  several  levels  of 
discretization  N\  lower  that  M(e),  for  which  Ajv,  is  also  smaller  than  e.  Let 
us  denote  by  N[{e)  the  lowest  level  of  refinement  satisfying  ^N[{e)  ^  The 
above  result  means  that  the  small  scales  can  be  integrated  with  a  larger  time 
step,  even  if  their  time  scales  are  smaller  than  the  ones  of  the  large  eddies. 
By  definition  and  due  to  the  behavior  of  |  Zjv,  I2,  with  respect  to  Ni,  the 
quantity  A/v,  increases  for  decreasing  values  of  N-i.  Hence,  the  appropriate 
time  step  for  iNi(e)  is  larger  than  the  appropriate  one  for  (n{{c)-  Moreover, 
on  Figure  4,  we  note  that  |  Zyv,  I2  presents  strong  variations  during  the  time 
evolution.  Consequently,  the  levels  Ni{e)  and  N[{e)  are  not  constant  in  time. 
To  take  into  account  this  specific  dynamics  of  the  small  scales,  we  propose, 
as  it  was  done  in  [10],  [7]  and  [17],  to  use  multigrid  technics.  The  integration 
of  the  intermediate  scales  ,  ^A/,(e)  <  ^Ni  <  consists  in  performing 

a  succession  of  multigrid  V-cycles  during  which  some  parts  of  the  spectrum 
are  frozen  during  the  time  evolution.  Hence,  we  use  here  a  property  of  the 
first  order  term  of  the  evolution  equation  of  the  small  scales. 

On  Figures  11,  we  can  see  that  the  variations  of  PjVi^mt(yAri>ZAf,)  and  of 
ZN^  are  correlated.  Therefore,  we  deduce  that  the  variations  of  Pv,  Pmt(y/Vi  >  ^Ni ) 
over  one  time  iteration  is  smaller  than  PutBintiVNu^Ni)  ?  fact,  one  can 
show  that  : 

Ajv,  =  I  PiVi  Pint(yyv, ,  Zw, )  I2  <  Cio  I  PyViPmt(y//,,ZAr,)  I2  • 

The  nonlinear  interaction  term  P/.  Pint(yNi  >  )  represents  a  small  part  of 

the  large  scales  time  derivative  .  Hence,  its  variation  over  one  time  itera¬ 
tion  can  be  neglected  without  a  lost  of  the  accuracy  on  the  large  scale  approx¬ 
imation.  We  use  here  a  property  of  the  time  evolution  of  Pv,  Pint(yAri  >  Zyv, ) 
which  is  a  second  order  term  of  the  large  scale  evolution.  As  the  time  scale  of 
ZjVj  is  much  smaller  than  the  one  of  ,  the  time  derivative  Pn,  Pint(yiVi ,  zyvi ) 
behaves  as  z^v,  and  then  is  very  oscilating  with  the  time.  From  the  previous 
inequality,  it  appears  that  the  quantity  A^,  can  be  controlled  by  the  term 
Pat,  Pint (y^Vi ,  zjv, )  ;  hence,  we  estimate  the  level  N[{e)  by  imposing  that  the 
ratio 

I  PjV>(e)Pint(yN((e)>ZA^>(e))  [2 
I  Pjv'(e)P(yN;(e),yyv;(e))  I2 
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remains  smaller  than  a  given  constant.  As  the  computation  of  the  nonlinear 
terms  requires  several  calls  to  the  FFT  routines,  the  evaluation  of  the  above 
ratio  is  expensive  and  then  prohibitive.  So,  we  have  looked  for  another  form 
of  it.  In  the  numerical  results,  we  have  noted  that 

I  PN^Biat{yN^^^N^)  U  ^  IIztVi ||  _  / ens(zAri) 

I  ,  y/v, )  I2  llyiv,  II  V  ens(yA/, ) ) 

as  we  can  see  on  Figures  13.  Hence,  we  use  the  ratio  of  the  enstrophy  of  the 
small  scales  over  the  enstrophy  of  the  large  ones  to  evaluate  the  level  N[{£). 
As  we  have  seen  before,  the  notion  of  variation  of  the  scales  is  intrinsic  to  the 
definition  of  N{{£).  More  sophisticated  criteria  will  be  derived  in  Section  2.2. 
Finally,  we  want  to  note  that  controlling  the  size  of  Ps^  Bii,t(yjv, ,  zn,  )  and  its 
time  variation  ensures  that  the  interaction  of  the  small  scales  over  the  large 
ones  is  negligible,  in  the  sense  that  this  interaction  can  be  locally  neglected 
from  a  numerical  point  of  view. 

2  Description  of  the  multiscale  method. 

As  it  was  shown  in  Section  1.2.2,  the  small  scales  of  the  flow  as  well  as  the 
terms  involving  their  nonlinear  interactions  can  be  fixed  in  time  during  a  few 
time  steps.  Nevertheless,  their  order  of  magnitude  may  change  draistically 
over  a  period  of  time ;  so,  the  cut-off  value  Ni  defining  the  separation  between 
the  small  and  the  large  eddies  can  not  be  let  fixed  in  time.  Hence,  we 
propose  a  multilevel  adaptative  procedure  evaluating  the  appropriate  level 
of  refinement  as  time  evolves.  Using  partly  theoretical  arguments,  we  show  in 
this  section  that  we  can  derive  estimates  for  the  variations  of  the  small  scales 
and  of  the  transfer  terms  to  the  largest  scales.  We  then  deduce  estimates 
for  the  length  of  the  frozen  periods.  Finally,  we  introduce  these  estimates  in 
the  algorithm  and  we  derive  a  dynamic  procedure  allowing  an  a  priori  and 
an  a  posteriori  control  of  the  length  of  the  quasi-static  time  intervals.  In  the 
first  subsection,  we  describe  the  multilevel  adaptive  procedure  ;  secondly,  we 
derive  time  scales  estimates  of  the  fine  structures  and  of  the  transfer  terms. 
In  the  third  part,  we  give  an  explicit  approximation  law  used  to  compute 
the  small  scale  components  of  the  flow  and  we  derive  error  estimates  for  this 
Approximate  Inertial  Manifold.  Finally,  we  describe  the  whole  algorithm 
including  the  modifications  previously  discussed. 
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2.1  The  multilevel  adaptive  procedure 

As  in  the  preceeding  Section  1.2.1,  we  choose  an  integer  N  larger  than  mo 
(cf  (8)),  which  represents  the  total  number  of  modes  retained  in  the  trunca¬ 
tion  and  we  denote  by 

k€/w 

an  approximation  of  u.  We  associate  with  N  the  largest  wave  number  of 
Uff,  namely  ks  =  N/2.  Hence,  the  smallest  scale  in  the  computation  is 
is  =  l/fc/v. 

We  are  now  given  a  sequence  of  levels  Ni  satisfying 

Ni  <  N2  <  ...  <  Ni  <  Afj+i  <  ...<  N.  (19) 

As  we  want  to  perform  pseudo-spectral  approximations  of  equations  (12)  on 
these  different  levels  of  refinement,  the  elements  Ni  of  these  sequences  have 
to  satisfy  the  restrictions  imposed  by  the  Fast  Fourier  Transforms  (FFT)  ; 
namely,  the  A^j’s  must  be  of  the  form  2’’  3’  S’",  where  p  >  2  and  9,  r  >  0.  Such 
an  algorithm  enables  us  to  define  a  suitable  sequence  of  levels  Ni.  Examples 
of  sequences  will  be  given  in  Section  3,  where  the  numerical  results  will  be 
described.  We  note  that  an  FFT  allowing  only  decompositions  in  powers  of 
2  is  not  efficient  for  this  purpose. 


Let  us  now  assume  that  the  approximation  u;v(x,f)  is  known  at  a  time 
tj.  As  it  was  suggested  in  the  previous  section,  we  define  two  levels  of  dis¬ 
cretization  Ni^{tj)  and  Ni^{tj)  by  the  following  procedure  : 


ft'i  is  defined  by  the  condition  that 
for  every  t  >  ti, 


ens(yN,ltj)) 


<0i. 


t2  is  defined  by  the  condition  that 
for  every  t  >  tj, 


(20) 


(21) 
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where  9i  and  are  two  given  constants  ;  62  is  chosen  so  that  e{zNij)  of 
order  of  the  accuracy  e*.  In  the  previous  versions  of  the  algorithm,  the  pa¬ 
rameter  $1  was  arbitrary  fixed  ;  we  will  derive  an  estimate  of  61  in  section  2.4. 
.  In  order  to  be  sure  that  Ni^  <  Ni^ ,  we  impose  the  additional  condition  that 
(*2  ~  *1)  bas  to  be  larger  than  a  given  constant.  This  is  motivated  by  the  fact 
that,  as  we  described  in  the  previous  section,  the  intermediate  region  of  the 
spectrum  between  and  is  a  transition  zone  between  the  large 

scales  and  the  static  (frozen)  small  scales. 

We  now  introduce  the  quantities 


Let  Tc{tj)  be  the  length  of  the  time  interval  during  which  the  scales  smaller 
than  ,  i.e.  ,  can  be  frozen  without  loosing  the  order  of  approximation 
on  the  larger  scales  ;  estimates  for  the  characteristic  length  Tc{tj)  will  be 
derived  in  Section  2.2.  In  previous  works,  see  for  instance  [8],  Te{tj)  was 
estimated  by  the  characteristic  relaxation  time  of  the  viscous  term,  namely  : 

Tc{tj)  = 

As  we  will  see  in  Section  3,  this  estimate  is  not  fine  enough  and  may  induce 
strong  errors  in  the  approximation  of  the  velocity  field.  The  available  levels 
of  refinement,  lying  between  and  are 

Ni^  <  Ni^+l  <...<Ni<...<  <  Nij, 

which  correspond  to  (*2  —  t'l  -1- 1)  levels.  For  the  sake  of  simplicity,  we  omit 
here  the  dependence  on  for  the  levels  Ni^  and  Ni,. 

As  in  classical  multigrid  methods,  we  use  the  concept  of  V-cycle  to  per¬ 
form  the  integration  of  (12)  on  the  interval  +  Let  us  define  a  V-cycle 
starting  at  time  tj.  Such  a  V-cycle  is  constituted  of  two  phases  described  as 
follows  : 
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•  phase  1  :  on  the  interval  +  (tj  —  zi)At],  the  current  level  Ni{t)  is 
defined  by  : 

Ni{t)  =  for  j  =  0, t2  -  ti, 

hence  Ni{t)  decreases  from  to 


•  phase  2  ;  on  the  interval  [tj  +  (tj  —  ii)At,tj  +  (2(t2  —  *i)  +  l)At],  the 
current  level  Ni{i)  is  defined  by  ; 

Ni{t)  =  for  j  =  ^2  -*i  +  I,...,2(i2-M)  +  1, 

hence  Ni{t)  increases  from  to 


Level  N, 


Then,  a  V-cycle  consists  in  [2(t2  —  ii)  + 1]  time  iterations.  The  quantity 
is  adjusted  so  that  the  time  interval  (t„  tj  +  Tc]  can  be  divided  into  a  fixed 
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number  of  such  cycles.  Figure  2  summarizes  this  process. 


Let  t  be  an  intermediate  time  on  the  interval  \tj,  tj  -{-  Tc]  ;  au:cording  to 
the  previous  procedure,  the  current  level  Ni{t)  is  given  by  : 


f  if  1  <r  <{t2-ti), 


mt)  =  { 


Ni. 


(24) 


if  (ij -*i  +  1)  <  r  <  2(t2  -  ii), 

where  r  is  given  by  : 

t-tj  =  (2p  (*2  -  *i)  +  r)  At. 

Knowing  the  size  Ni  of  the  coarse  grid  at  time  t,  we  decompose  Ujv(t)  into  : 

uyv(t)  =  yNi{t)  +  ZNAt), 

where  ysiit)  represents  the  scales  larger  than  and  the  scales 

smaller  than  tsi  and  larger  than  is-  The  computation  of  both  components 
ySiii)  and  ZSi(t)  are  performed  as  follows  : 

•  computation  of  zsi(t) : 

ZAr^(t)  =  zsi(t-At), 
i.e.  zsi{t)  is  frozen  and  set  to  its  last  value. 

•  computation  of  ysM  '  in  order  to  evaluate  ysi(i),  we  integrate 
equation  (12)  over  the  interval  {<  —  At,  <]  ;  then  it  follows  that 

ti(k,t)  =  c-‘'W*^‘ti(k,t-At) 

/-A*  yyVi(T))  dr  (25) 

+  /  (t),  m ('t))  dr 

for  every  k  in  /jv,  =  [1  —  Nif^yNij2]  x  [0,  iV,/2]. 
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The  first  integral  is  computed  by  an  explicit  Runge-Kutta  scheme  of  order 
3  (see  [7]).  With  this  scheme  the  interval  [t  —  At,t]  is  split  into  three  sub¬ 
intervals  of  the  form  [tj,  <i+i],  where  to  =  t  —  At  and  ^3  =  t.  On  each  of  these 
sub-intervals,  the  second  integral  is  approximated  as  follows  ; 

•  (26) 


We  note  here  that  this  approximation  is  an  explicit  Euler  scheme  on  the 
interval  where  the  following  approximation  is  performed 


This  integration  requires  the  storage  of  PN,BM{yNi{ti),^N,{tj))f  at  the  be¬ 
ginning  of  the  cycle,  for  each  coarse  grid  Ni  between  Ni^  and  Ni^,  i.e.  for 
(*2  —  *1  +  1)  levels. 


As  the  time  Te{tj)  is  adjusted  to  be  a  multiple  of  a  complete  V-cycle, 
the  current  level  at  time  tj  -f-  Te{tj)  is  equal  to  Aj,,  the  highest  coarse  level. 
Referring  to  the  integration  process  described  above,  the  large  scales  y^^^ 
are  known  at  the  end  of  the  cycle,  i.e.  at  time  tj  -f  Tc{tj).  The  smallest  scales 
are  then  updated  by  projecting  the  approximate  solution  u/v(<j  -l-Tc(tj)) 
on  an  approximate  form  of  the  small  scales  equation  (12),  for  instance 

+  Tc{tj))  =  <i>iyNi^{ti+rcitj)),ZNi^{tj),Tc{tj)).  (27) 

(27)  is  the  equation  of  an  Approximate  Inertial  Manifold.  Such,  manifolds 
were  first  derived  in  [1]  ;  (27)  provides  an  interaction  law  between  the  small 
and  the  large  scale  components  of  the  flow,  and  expresses  as  a  function 
of  yNi2‘  Further  information  on  such  law  can  be  found  in  [isf,  [19],  and  [7]. 
Other  kinds  of  approximate  inertial  manifolds  are  derived  in  [20]  and  [21]  ; 
the  implementation  of  an  Approximate  Inertial  Manifolds  of  first  order  will 
be  discussed  in  Section  2.3.  In  Section  2.3,  we  will  discuss  on  the  eflSciency 
of  these  nonlinear  forms  and  derive  some  estimates  which  explain  in  which 
range  of  the  spectrum  they  can  be  implemented.  From  a  strictly  compu¬ 
tational  point  of  view,  an  approximate  inertial  manifold  is  efficient  in  the 
sense  that  this  equation  allows  us  to  estimate  the  small  scales  as  an  explicit 
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function  of  the  large  ones.  Moreover,  a  first  order  law  only  requires  one  eval¬ 
uation  of  the  nonlinear  terms,  on  the  fine  grid  (see  Section  2.3). 

Finally,  at  tj+j  =  tj  +  Tc(<j),  i.e.  at  the  end  of  a  whole  cycle,  we  have  an 
approximation  of  UN{tj+\)  by  : 

un(<j+i)  =  y;v.j(<j+i)  + 

We  start  the  full  procedure  again  by  computing  two  new  levels  and 

Then,  we  perform  new  V-cycles  on  the  time  interval  [tj+\,tj+i  + 
Tc(tj+i)],  and  so  on. 


Basically,  we  can  summarize  this  process  by  saying  that  the  full  time  in¬ 
terval  [0,  T]  of  the  whole  computation  is  split  into  several  small  time  intervals 
-1-  Tc{tj)\  and,  on  each  of  these  intervals  the  velocity  spectrum  is  split 
into  three  fundamental  regions  : 

•  the  dynamical  range,  corresponding  to  wavenumbers  smaller  than  . 
It  represents  the  large  scale  structures  of  the  flow  -  i.e.  the  scales  con¬ 
taining  most  of  the  energy  and  the  enstrophy  of  the  flow.  The  modes 
lying  in  this  part  of  the  spectrum  are  integrated  at  each  time  step  of 
the  time  interval  +  Tc{tj)]. 

•  a  transition  range,  corresponding  to  wavenumbers  between  kN.^  and 

An  up  and  down  oscillation  process  is  used,  i.e.  the  current  level 
of  discretization  Ni{t)  undergoes  all  the  intermediate  levels  between 
and  while  the  time  evolves. 

•  a  quasi-static  range,  i.e.  for  wavenumbers  larger  than  k^-^.  It  repre¬ 
sents  the  smallest  scales,  which  are  numerically  negligible  ;  i.e.  their 
energy  and  their  variations  are  smaller  than  the  expected  accuracy. 

Figure  3  represents  these  three  different  regions  of  the  spectrum. 

Now  we  want  to  make  some  remarks  on  the  effect  of  the  integration 
procedure  previously  described  on  the  smale  scale  components  of  the  flow. 
One  of  the  crucial  points  concerns  the  technic  used  to  update  the  small 
structures  Zjv,  (0  lying  in  the  transition  range. 
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Dynamic 


E(k) 


Figure  3:  The  multilevel  procedure. 

Let  us  first  recall  the  equation  satisfied  by  the  k—  coefficient  of  the  Fourier 
decomposition  of  uyy  : 

+  1/ 1  k  p  u(k)  +5k(yAr, ,  yn,  ) 

+  5int,k(yNnZN,)  =  0. 

We  assume  here  that  k  is  such  that  |  k  |>  mo  and  hence  (see  (8))  : 

g(k)  =  0. 

We  can  then  write  (28)  in  the  following  form 


(28) 


—  (e‘'*’‘u(k))  =  .  (29) 


Integration  of  (29)  over  a  time  interval  [t,  t  +  r]  leads  to 
u(k,  t  +  T)=  e“‘'l*'**’’u(k,  t) 

t+T 


(30) 


We  note  that  if  k  is  large  enough,  u(k)  is  a  component  of  the  small  scales. 
For  instance,  we  assume  that  : 

Then,  k  lies  either  in  the  transition  or  quasi-static  range.  According  to  the 
implementation  of  the  multilevel  method  in  (30),  u(k,t)  is  replaced  by  an 
approximation  close  to  the  actual  value,  as  the  perturbations  are  smaller 
than  the  accuracy  £.  If  k  lies  far  enough  in  the  dissipation  range,  then  the 
dissipative  terms  are  quasi-dominent,  namely  : 

I  Z/v,  I2  ~  1  l/AZNt  I2  I  Qjv,^mt(yjv,,zjv,)  I2, 

as  we  can  see  on  Figure  12.  In  that  case,  the  errors  introduced  by  the  quasi¬ 
static  integration  are  quickly  damped  by  the  effects  of  the  operator  —A  : 
few  corrected  iterations  are  needed.  If  k  is  not  as  large,  i.e.  if  it  is  closer  to 
the  inertial  range,  then  the  coupled  nonlinear  terms  become  the  most  impor¬ 
tant  terms  of  equation  (30).  In  that  case,  the  error  is  convected  by  nonlinear 
effects  in  the  largest  wavenumbers,  and  are  finally  damped  by  viscous  effects. 

During  the  V-cycles,  the  modes  closer  to  iV„  (t)  are  more  often  integrated 
than  the  ones  close  to  Ni,{t).  The  structure  of  the  V-cycles  is  thus  well- 
adapted  to  the  integration  of  the  small  scales.  As  we  will  see  in  Section  3, 
devoted  to  the  description  of  the  numerical  results,  there  is  no  accumulation 
of  errors  in  the  intermediate  scales  and  there  is  no  energy  pile-up  in  the  high 
wavenumbers  of  the  Fourier  decomposition.  The  enstrophy  cascades  axe  well 
described  by  the  V-cycle  technic. 


2.2  Time  scale  estimates  for  the  small  eddies  and  for 
the  nonlinear  interaction  terms. 

Time  scale  estimates  for  Zjv, . 

From  now  on  and  for  the  sake  of  simplicity,  we  omit  the  subscripts  Ni  in  the 
notations.  We  then  denote  by  z{t)  the  small  scale  components  of  the  flow. 
As  the  Navier-Stokes  equations  are  analytic  in  time  (see  [1]),  we  can  write  a 
Taylor  expansion  of  z{t)  : 

z{t  +  T)  =  z{t) TZ{t)  +  o{t^).  (31) 
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As  we  saw  in  Section  1.2.2,  the  quantity  rz(t)  can  be  much  smaller  than  z(t). 
We  can  then  assume  that  the  higher  order  terms,  in  the  Taylor  expansion 
are  negligible  by  comparison  with  the  first  order  ones. 

As  in  the  previous  section,  we  denote  by  e  the  accuracy,  i.e.  the  solution 
is  approximated  with  an  error  of  the  order  of  e.  Let  us  assume  that,  on  an 
interval  {<,  <  +  t],  we  tolerate  an  error  of  the  order  of  e  on  the  approximation 
of  the  small  scale  components.  From  (31),  we  then  derive  an  estimate  on  t, 
namely  : 

^  Ke 

^-\m  b' 

We  denote  t(£)  =  Ke/\  z(t)  jj,  where  is  a  nondimensional  real  constant  of 
the  order  of  the  unity.  As  the  time  derivative  |  z(t)  I2  has  a  decaying  spec¬ 
trum,  the  quantity  T(e)  then  decreases  with  decreasing  values  of  the  level 
Ni.  Estimate  (32)  then  provides  a  restriction  on  the  available  level  on  which 
the  modes  can  be  frozen,  i.e.  there  exists  a  level  Ni  such  that  t(£)  becomes 
smaller  than  the  time  step  At.  On  Figure  14,  we  have  plotted  the  ratio  r(e) 
for  different  values  of  Ni.  These  results  are  obtained  from  the  computation 
presented  in  Section  3.  In  this  case,  the  accuracy  is  given  by  the  temporal 
discretization  and  is  of  the  order  of  At^.  As  the  time  step  At  is  of  the  order 
of  10~^,  we  can  estimate  Ke  ~  10“*.  So,  for  a  fixed  given  value  of  the  level 
Ni,  the  corresponding  time  scale  T{e)  presents  very  strong  variations.  In  a 
previous  version  of  the  algorithm,  t  was  estimated  by  and  then 

was  constant  in  time  ;  this  choice  is  obviously  inappropriate.  Deriving  an 
estimate  on  r(e)  is  essential  to  insure  the  efficiency  of  the  algorithm.  Indeed, 
we  can  imagine  a  situation  where  the  procedure  (20)  and  (21)  allows  the 
choice  of  a  level  TV,-,  while  the  condition  (32)  is  violated,  i.e.  that  r  is  much 
smaller  than  the  time  step.  Hence,  the  constraint  (32)  provides  an  additional 
way  to  determine  Ni^  and  Ni^. 

In  the  transition  range  [A^»,(<j),  A/i2(tj)],  we  have  a  time  estimate  TjVi(e), 
given  by  (32),  for  level  Ni  : 


'TNiie) 


Ke 


As  we  have  noted  above,  the  quantity  r^tie)  decreases  when  N  decreases. 
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According  to  the  definition  of  a  V-cycle,  on  the  lower  level  the  char¬ 

acteristic  time  scale  h2LS  to  satisfy  : 

TN,^ie)>At.  (33) 

We  recall  that  (33)  is  motivated  by  the  fact  that  the  scales  will  be 
frozen  over  only  one  time  iteration  during  a  complete  V-cycle  ;  (33)  is  then 
a  constraint  that  has  to  satisfy.  As  we  have  previously  remarked,  the 

time  derivative  can  be  evaluated  by  the  nonlinear  terms,  i.e.  : 

la  I  QNi,B{UN,Us)  U  ■ 


At  time  tj,  B,n«(yjv.,  >  *JVi, )  is  obtained  by  using  the  following  relation  : 
fiint(yyv,,,ZN,,)  =  B{un,un)  - 

Hence,  B{un,  ua^)  will  be  computed  at  that  time.  Then,  the  computation  of 
the  quantity 


Ke  ^  Ke 


(34) 


will  not  add  extra  cost.  With  (33),  we  are  sure  that  on  all  levels  iV,  higher 
than  Ni^ ,  the  corresponding  scales  can  be  frozen  during  more  than  one  time 
iteration. 


The  characteristic  time  TNi^{e)  provides  an  estimate  of  the  global  time 
length  Tc{tj)  of  the  whole  cycle  -j-  Tc{tj)]  ;  7Wi,(e)  can  be  evaluated  as 
in  (34),  i.e.  : 

I  1^- 

Remark  1  :  if  Ni^  lies  in  the  quasi-static  range,  another  estimate  can  be 
derived  ;  as  we  have  seen  before,  we  have  : 


I  I' \  I2  . 

Then,  we  can  write  : 


1  zjVi,(^j)  I2  >  cnv  I  AzNi^{tj)  I2, 
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where  Cn  is  a  constant  of  the  order  of  unity.  FVom  the  definition  of  the  norm 
I  —  I2  and  of  zjVij ,  we  can  obtain  ; 

1/  1  AZNi^itj)  \i>  V  I  ZN^{tj)  I2  . 

We  then  obtain  the  following  estimate  of  t^v*,  (e)  : 

(  .  _ _ 

I  Ij' 

Moreover,  we  recall  that,  by  definition  : 


We  finally  have  an  a  priori  estimate  for  ta/.^  : 


_ Ke _ 

1/  kN^%{tj)  1  I2  ■ 


(36) 


If  £  is  the  accuracy  of  the  time  scheme,  we  recall  that  K  corresponds  to  a 
high  order  derivative  of  U;v  versus  time.  Hence,  we  can  reasonably  assume 
that  K  is  at  least  of  the  order  of  |  yNi^{tj)  I2  •  This  case  occurs  when  a  Direct 
Numerical  Simulation  is  performed. 


Time  scale  estimate  for  the  transfer  terms. 

Let  us  denote  by  y  instead  of  y^^.  the  large  scale  component  of  the  flow.  We 
recall  that  y  is  governed  by  the  equation  : 

y-t/Ay+ PB(y,y)  + PBint{y,z)  =  Pg. 

Here,  P  denotes  Pff^  and  y  =  Let  us  rewrite  the  previous  equation  under 
the  form  : 

y  =  uAy  +  Pg-  PB{y,y)  -  PBint{y,z). 

We  introduce  the  function  F  defined  by 

F(y)  =  uAy  +  Pg-  PB{y,y), 

so  that  : 

y  =  F(y)-PBi„t(y,z).  (37) 
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In  this  form,  it  clearly  appears  that  the  transfer  terms  PBintiViZ)  is  a  cor¬ 
rection  to  the  time  derivative  of  the  large  scale  components.  As  it  was  done 
before  for  the  small  structures,  we  derive  a  Taylor  expansion  of  y(<)  : 

y(<  +  T)  =  y(0  + ’■y(0  +  yy(0  + ‘Kt®).  (38) 

We  assume  here  that  the  terms  of  order  larger  than  three  are  negligible. 
From  (37),  we  can  derive  : 

y(0  =  ^(y)-Mnt(y,z). 


Reporting  this  expression  into  (38),  we  obtain  : 

y(t4T)  =  y(0  +  Ty(0  + Y^(y)- Y^^.nt(y,z)  +  o(T3). 


As  we  tolerate  an  error  of  the  order  of  e  on  y{t  +  t),  the  coupled  nonlinear 
terms  PBint(y,  z)  can  be  frozen  during  a  time  t,  if  : 

j  \  PBin, iy,z)\2  <  Ke.  (39) 

Then,  it  follows  an  estimate  on  r  : 


/  2Ke 
\|  Ffimi(y,z)  I2  j 


r'(e). 


(40) 


We  note  that,  if  e  is  the  accuracy  of  the  scheme,  condition  (40)  is  necessary 
to  preserve  the  order  of  the  time  scheme.  On  Figure  15,  we  have  plotted 
the  evolution  of  the  ratio  T'{e)  for  different  levels  of  refinement  AT,-.  As  for 
T(e),  the  quantity  T'{e)  decreases  when  Ni  decreases,  which  is  due  to  the 
fact  that  PBint{y,z)  has  a  decaying  spectrum,  like  z(<).  So,  for  the  levels 
Ni  lying  in  the  transition  range,  i.e.  between  Ni^  and  Ni^,  the  vadue  of  T'(e) 
corresponding  to  the  level  ATjj  is  the  most  restricted  one.  Hence,  in  order  to 
control  the  variations  of  PNiBiatiysn^Ni)  on  the  different  levels,  a  sufficient 
condition  is  to  estimate  P{e)  on  the  lower  level  A/,,  of  the  transition  range. 
We  want  to  note  that  the  mathematical  estimates  which  can  be  derived  on 
the  time  derivative  PBint{y,z)  of  the  transfer  terms  do  not  provide  efficient 
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information  ;  nevertheless,  the  numerical  experiments  show  a  correlation 
between 

1  pgmt(y.z)  u  ^  ,  I  m  b 

1  /'5int(y,z)  I2  |z(0|2 

as  we  can  see  on  Figure  11.  So,  we  deduce  that 

y  I  PBintiy,z)  I2  >  Ci2  Y  I  PBint(y,z)  I2 


where  c  is  nondimensional  constant  of  the  order  of  the  unity.  Hence,  it  follows 
that  : 


2Ke 


z{t) 


C12  I  Z(«)  I2  I  Pfiint(y,z)  I2 


r= 


t"(£). 


(41) 


We  then  obtain  with  (41)  an  estimate  of  t  as  a  function  of  ]  i{t)  I2  .  From 
a  computational  point  of  view,  (41)  is  much  more  eflficient  than  (40).  As  it 
was  noted  before,  we  have  to  derive  an  estimate  on  rjy.^  (e)  in  order  to  control 
the  time  variations  of  the  transfer  terms  which  depend  on  the  scales  in  the 
transition  range.  We  now  recall  that  the  level  A/j,  is  defined  by  the  evaluation 
of  the  ratio  : 


llyN..(<i)ir 


Moreover,  this  ratio  is  equivalent  to  the  ratio  of  the  coupled  nonlinear  terms 
of  the  large  scales  in  the  energy  norm,  namely 


I  PNi^  ^int(yAfi, ,  )  I2  ^  II  II 

I  Psi^  ,  y;v,, )  I2  ~  1 1  y^v.,  (< > )  1 1 


Cl3 


where  c  is  of  the  order  of  the  unity.  The  estimates  can  then  be  written  under 
the  new  form  : 


2Ke 


_  ^  I  {tj)  I2 

)  I  Pn,,  B(yN,, ,  ySi, )  I2  I  isi,  {tj)  I2 


) 


1/2 


(42) 


In  (42),  the  time  derivative  of  the  small  scale  components  can  be 

estimated  as  it  was  done  previously.  Finally,  TSi^  (e)  provides  a  constraint  on 
the  choice  of  the  level  TVi,,  while  Tff^{€)  and  provides  two  estimates 

of  the  length  of  the  whole  cycle. 
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2.3  Approximate  equation  for  the  quasi-static  scales 
:  projection  on  an  Approximate  Inertial  Man¬ 
ifold. 

In  this  section,  we  intend  to  discuss  the  efficiency  of  the  Approximate  Inertial 
Manifolds  (AIM)  and  show  in  which  part  of  the  spectrum  they  can  be  used. 


We  consider  an  approximation  of  the  equation  of  the  small  scales  (12)  in 
which  we  drop  the  coupled  nonlinear  terms  : 

~-i^Az  +  QB(y,y)  =  0.  (43) 

We  now  introduce  an  operator  defined  by  : 

^  c*'W'‘ti(k,f)u;k- 

k€/w\/w, 

Then  by  applying  this  operator  to  (43),  we  obtain  : 

=  -e-«'^0B(y,y).  (44) 

We  assume,  in  agreement  with  the  method  described  in  Section  2.1,  that  z 
has  been  frozen  on  the  time  interval  +  Tc{tj)]  where  Tc{tj)  is  estimated 
as  in  the  previous  section.  For  the  sake  of  simplicity,  we  write  here  t  instead 
of  tj  and  Tc  instead  of  Tc(tj).  We  then  integrate  (44)  over  the  interval  [f,  f+rj, 
which  yields  : 


z(f  +  Tc)  =  e*'"‘^z(<)- J^‘'^^‘e-'(‘'-‘-"*)^QB(y(<r),y(cr))d(T. 
Consider  then  the  following  approximation  of  the  right-hand  side  : 
e-"<‘'-(‘+^‘»^gB(y((r),y(<r))  da 

-  -  e'^‘^)Q5(y(<  +  rc),y(< +  Tc)). 

Finally,  z{t  +  Tc)  is  computed  by  : 
z(<  +  Tc)  =  c*^‘=^z(t) 

-  e'"‘^)QB{y(t  +  T,),y{t  +  rj) 


(45) 


(46) 
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(46)  means  that  the  smal]  scales  are  slaved  by  the  large  ones.  Also,  from  a 
computational  point  of  view,  (46)  provides  an  efficient  way  to  evaluate  the 
small  scales.  In  comparison  to  a  classical  time  scheme,  (46)  presents  the 
advant2ige  that  only  one  evaluation  of  the  nonlinear  terms  is  required.  The 
error  occuring  by  using  (46)  instead  of  integrating  the  small  scales  equation 
is  constituted  by  two  components,  namely  the  time  discretization  on  the 
approximation  of  the  integral  and  the  dropped  terms  QB^niy,  z)  : 

I  £(*)  ll  <  I  z(»))  da  |, 

+  I  l‘*'‘  (QB(y(<T),  y(a))  -  (?B(y((  +  r.),  y{t  +  r.)))  da  |, 

<  Tc  I  <3B«t(y,z)<.rj2  +  2  Te  1  QB(y ,  y)<.Tj2 

(47) 

where  |  QBintiy,  x)t,rc  I2  =  max  |  QBiat(y(<y),  z{<t))  I2,  and  similarly  for 

<y6l<,«+Te] 

1  QB{y,y)t,rc  I2 . 

Considering  the  previous  discussions, 

|a:(<)|2  ~  lQBi„t(y(f),z(<))l2 

|<?^int(y,z)t,Tj2 

^  1  Q5(y,y)t.rc  I2 . 

Then,  we  obtain  an  estimate  of  the  error  : 

1  £(z)  |2  <  Ci4  Tc  I  z(<)  I2  .  (48) 

Recalling  that  Tc  is  defined  such  that 

Tel  i{t)  I2  <  K  e, 

we  then  have  the  following  estimate  : 

I  e(z)  I2  <  cj4  K  £. 

Then,  from  the  definition  of  the  level  Ni^  and  Tc,  the  error  introduced  by 
using  the  approximate  equation  (46)  is  always  smaller  than  s.  We  want  to 
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note  that  the  error  e(z)  does  not  depend  directly  on  the  level  Ni^  :  there  is 
no  restriction  on  the  value  of  iVj,.  The  Approximate  Inertial  Manifold  (46) 
can  then  be  used  to  simulate  the  evolution  of  the  fine  structures  of  the  flow 
even  for  wave-numbers  lying  inside  the  inertial  range  of  the  spectrum.  In 
the  numerical  simulations  presented  in  the  next  subsection,  (46)  will  be  used. 

Remark  2  :  Let  us  now  consider  the  Approximate  Inertial  Manifold  intro¬ 
duced  in  [1],  namely  : 


-  t/Az  +  QB{y,y)  =  0.  (49) 

We  recall  that  (49)  consists  of  an  approximation  of  the  ftill  equation  of  the 
small  scale  components  z,  where  the  time  derivative  z  as  well  as  the  coupled 
nonlinear  terms  QBint(y,z)  have  been  dropped.  In  the  case  considered  here 
of  periodic  boundary  conditions,  it  is  easy  to  invert  the  Stokes  operator  (—A) 
and  then  z  can  be  evaluated  as  a  function  of  y  : 

z  =  {-vA)-'QB(y,y).  (50) 

At  this  point,  we  note  that  for  large  values  of  Tc,  (46)  and  (50)  are  equivalent. 
The  order  of  magnitude  of  the  dropped  terms  in  the  small  scde  equations 
may  induce  a  restriction  on  the  use  of  (50)  to  evaluate  z.  In  fact,  we  want 
to  find  criteria  telling  in  which  range  of  the  spectrum  (50)  can  be  applied. 
The  spatial  error  ^(z)  appearing  when  (50)  is  used  is  exactly  given  by  the 
difference  between  (50)  and  the  z  equation,  i.e. : 

5(z)  =  (i/A)“*(z-|-QBint(y,z)). 

We  can  then  estimate 

I  5(z)  la  <  I  i  +  QBmt(y,z)  la  • 

In  the  numerical  experiments  that  we  have  conducted,  we  have  seen  that  the 
right-hand  side  of  the  previous  inequality  is  of  the  Order  of  (*^fc^,)~^  I  ®  |a  i 
hence  it  follows  that 


|^(z)|a<  cn{uk%j-^\z{t)\2, 

where  Cis  is  a  nondimensional  constant  of  the  order  of  unity.  Then,  the  spatial 
error  introduced  by  (50)  mainly  depends  on  the  size  of  |  z{t)  {2  •  Assuming 
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that  Ni  lies  in  the  dissipation  range  and  that  |  z(t)  I2  can  be  estimated  by 
I  i/Az  I2,  we  obtain  : 

I  Six)  I2  <  c,6  I  z  I2, 


by  using 

|zyv.  12^^  iuk%J  Izjv,  I2. 

Hence,  in  that  specific  case,  the  error  ^(z)  is  of  the  order  of  z  itself.  Thus  (50) 
can  be  used  in  the  quasi-static  range  of  the  spectrum  where  z  is  of  the 
order  of  the  scheme  accuracy.  On  Figure  16,  we  can  see  that  the  quantity 
I  I2  becomes  larger  than  |  z  I2  itself  when  the  cut-off  value  Ni 
decreases.  Hence,  it  seems  that  if  (50)  is  an  efficient  way  to  compute  the 
very  fine  structures  of  the  flow  lying  far  inside  the  dissipation  range,  it  is  no 
longer  the  case  for  the  scales  of  the  order  of,  and  immediately  larger  than 
the  dissipation  scale  if,. 


2.4  Description  of  the  complete  multilevel  algorithm. 

In  this  section,  we  summarize  the  complete  multilevel  method  which  includes 
the  time  scales  derived  in  Section  1,2.2.  We  still  denote  by  e  the  accuracy 
of  the  computation  ;  we  recall  that  £  is  a  given  parameter  in  the  following 
algorithm. 

As  in  subsection  2.1,  we  choose  a  sequence  of  levels  Ni  such  that  : 

Nx  <  Ni  <  . . .  <  Ni  <  Ni+x  <  ...  <  N. 

The  whole  time  interval  of  the  simulation,  namely  [<o»to  +  T],  is  split  into 
several  intervals  of  the  form 


where  tj  =  <0  +  Tc(tfc).  Futhermore,  we  assume  that  the  final  time  to  +  T 

k=l 

satisfies 

to  +  T  tm, 

so  that  we  have  m  intervals.  Let  us  now  assume  that  the  approximate  solution 
u/v(x,t)  is  known  at  time  tj  with  j  <  m.  As  in  2.1,  we  compute  two  levels 
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of  refinemeDt  Ni^{tj)  and  according  to  the  procedures  (20)  and  (21). 

Moreover,  we  impose  that 


Ke 

I  u 


>  At. 


With  (35)  and  (42),  we  derive  an  a  priori  estimate  of  respectively  (e)  and 
(e)  ;  we  then  obtain  an  evaluation  of  the  length  Tc{tj)  of  the  j—  cycle  : 

Tc(<j)  =  min(Tiv.j£),  t;J.^(£)).  (51) 

We  note  here  that,  with  this  definition,  Tc(tj)  can  be  smaller  than  one  V- 
cycle,  i.e.  (2(i2  — *i)  +  l)At ;  in  such  case,  levels  and  Ni^  are  too  small  and 
need  to  be  reevaluated.  With  (51),  we  have  an  a  priori  estimate  of  the  global 
length  of  the  whole  cycle.  Finally,  we  have  computed  the  three  characteristic 
values  : 

As  in  2.1,  the  integration  is  performed  on  the  interval  +  ’■c(ti)]  by  a 
succession  of  V-cycles.  At  the  end  of  each  V-cycle,  i.e.  at  time  tj+pv  = 
tj  +  (2p(t2  —  ti)  +  l)Af,  we  derive  an  a  posteriori  estimate  of  the  quantities 
T1v<,  (^)  ^Ar<j(«)«  Hence,  we  take  into  account  the  evolution  of  the  scales 

lying  in  the  transition  range  of  the  spectrum  of  the  velocity  (see  Figure  3).  At 
this  time,  if  (t>+Tc(tj))— fj+pv  is  larger  than  one  full  V-cycle,  i.e.  2(^2— *i)+l 
time  iterations,  we  perform  another  V-cycle  after  reajdusting  the  value  of 
Tc{tj)  with  the  new  values  of  (£)  and  Now,  in  the  other  case,  i.e. 

if  {tj  +  Tc{tj))  —  tj+pv  is  smaller  than  one  V-cycle,  we  stop  the  whole  cycle 
by  saying  that 

+  'rc(fi)  =  ^i+i* 

At  tj  +  'rc{tj),  we  compute  the  small  scale  components  of  the  spectrum 
by  projecting  the  solution  on  the  Approximate  Inertial  Manifold  (46). 

Then,  we  readjust  the  two  levels  and  restart  a  new  cycle  as  it  was  done  in 
Section  2.1. 

Before  we  conclude  this  section,  we  want  to  note  that  with  this  algorithm, 
in  opposition  with  the  previous  ones  (see  for  instance  [10],  [7]  or  [17]),  the 
constants  0i  and  62  of  the  procedures  (20)  and  (21)  can  be  fixed  very  easily. 


35 


Indeed,  as  we  have  previously  seen  the  parameter  is  chosen  so  that  c(za^,^  ) 
is  of  the  order  of  e*.  The  constant  di,  which  provides  an  estimate  of  the  ratio 
II  II  /  II  yjVi,  IK  can  be  evaluated  at  the  initid  state,  i.e.  t  =  0,  by 

“  Atfjy/v,,  |||y^,,  W 


Hence,  we  insure  that  : 

Atv,,  =  At  1  |2<  cAt  I  I2II  ||<  e. 

Then,  the  condition  on  tw,j(,)  is  satisfied  i.e.  ^  At.  Moreover,  we 

can  implement  a  self-adaptative  procedure  allowing  a  dynamical  reevaluation 
of  these  constants  during  the  time  evolution.  So,  if  0i  and  02  were  previously 
fixed  in  an  empirical  way  in  the  algorithms,  we  have  found  now  a  more 
efficient  way  to  evaluate  these  constants. 
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Figure  4:  Time  evolution  of  the  ratio  \ — for  Ni  =  32,  64,  128  and 

1  y^i  u 

M  =  196. 
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Figure  5:  Time  evolution  of  the  ratio  i  ^  for  N-i  = 

*  (  PN,B{yN,,yN,)  b 

32,  64,  128  and  Nx  =  196. 
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Figure  6:  Time  evolution  of  j  i^v,  I2  and  v  |  Azat,  I2  for  N]  =  128  and 
Nx  =  196. 
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Figure  7:  Time  evolution  of  |  zjv,  I2  and  |  Q%^Bint{yNi,^Ni)  I2,  for  N\  =  32 
and  N\  =  64, 
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Figure  9:  Time  evolution  of  |  isi  ia  for  =  32,  64,  128  and  A^i  =  196. 
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Figure  10:  Time  evolution  of  [  la,  v  |  AyN^  la,  |  ^(yw, ,  Xn,  )  |a  for 

yVi  =  32. 
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Figtire  11:  Time  evolution  of 
32,  64,  128,  and  196. 
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Numerical  results 


3.1  Comparison  with  the  Galerkiu  method  and  with 
a  previous  version  of  the  multilevel  method. 

In  this  section,  we  report  numerical  results  obtained  by  using  the  nonlinear 
Galerkin  method  and  a  classical  Galerkin  projection.  The  flow  of  Kolmogorov 
type  is  forced  by  a  time  independent  external  force  f  which  acts,  in  the 
spectral  space,  on  only  some  low  frequency  components  of  the  velocity  field. 
The  initial  condition  is  chosen  so  that  its  spectrum  has  a  specified  shape  but 
the  phases  of  its  Fourier  components  are  randomly  chosen.  So,  the  flow  at 
time  t  =  0  has  no  organized  struture.  We  have  let  the  flow  evolve  on  over 
10*  time  iterations,  i.e.  from  t  =  0  to  <  =  100  ;  this  is  much  longer  than  the 
integral  time  scale,  which  is  of  the  order  of  the  unity.  We  have  compared  the 
solutions  obtained  with  both  methods. 


3.1.1  Description  of  the  computation. 

The  initial  condition  here  is  computed  from  a  given  spectrum  of  the  initial 
vorticity  Uq  =  V  x  where  is  the  following  expansion  : 


=  S  Uo,k(<)e’^> 

ke/w 


with  X  =  (a:tia?2)  €  fl.  We  choose  by  setting 

Wo,k  =  I  «'o,k  I 

where  0^  €  [0, 2ir]  is  generated  by  a  random  function,  and  ; 

if  ib=lkl<fc„. 


(52) 


(53) 


‘*'o,k  I  — 


0 


(54) 


otherwise  ; 


Ci7  is  determined  such  that  1  ujq  2.0  ;  ka  is  equal  to  60.  At  this  point, 
we  note  that  if  (  u)o,k  H  energy  spectrum  of  Uo 

=  S  (I  «o,k  I*  +  I  Uo,k  I*) 

|k|  =  k 
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is  like  So,  in  the  case  presented  here  =  —0.5  and  we  have  Eoi^k)  ~ 

k~^.  On  Figure  28,  one  can  see  the  isovorticity  lines  of  the  initial  velocity 
field  and  Figure  17  shows  the  energy  spectrum  E{k). 

The  external  force  f  is  constant  in  time  and  has  only  a  few  non-zero 
wavenumbers,  namely  : 

4  =  (/i.k,/2.k), 

with 

I  l/i*l=  c,8  if  k€ZVl*i  l+l*:2l=  3, 

I  I  /i,k  I  =  0  otherwise, 

Ci8  is  determined  such  that  |  f  [2=  0.225.  The  Fourier  coefficients  of  f  are 
finally  obtained  by 

/a  =  I  /i*  I 

where  the  phases  €  [0,2ir]  are  randomly  generated. 

In  order  to  describe  all  the  scales  of  motion,  the  number  of  modes  N 
in  each  dimension  of  the  space  must  be  chosen  so  that  the  associated  grid 
size  2iflN  is  smaller  than  the  dissipative  (Kraichnan)  scales  i,,  ;  in  term 
of  wavenumber,  it  means  that  ks  >  K,.  We  recall  that  under  the  dissipative 
scales,  the  motion  is  damped  by  viscosity.  In  fact,  the  total  number  of  degrees 
of  freedom  needed  to  describe  the  motion,  from  the  dissipative  scales  to  the 
large  scales  containing  eddies,  can  be  estimated  by  the  ratio  {knlko)^  (see  [22] 
and  [16]).  Constantin,  Foias,  Manley  and  Temam,  in  [16],  have  related  this 
quantity  to  the  dimension  of  the  attrator  of  the  Navier-Stokes  equations 

where  Rei  is  the  integral  scale  Reynolds  number,  which  can  be  defined  by 


Here  i?*  =  (2/  |  Cl ))  c(u)  and  L  =  l/ki,\s  the  integral  length  scale  defined 
as  ; 

L  =  — where  E  =  (1/2)  and  7  =  i' )  Au  [2  /  j  j  {v  =  10“^). 
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L  is  of  the  order  of  1.72  and  d  is  of  the  order  of  0.45  (see  Figure  19).  It 
follows  that  Rei  =  784  and  kt,  =  17,  which  corresponds  to  N,,  =  2  kt,  =  M 
modes  in  each  direction  of  space.  In  order  to  be  sure  to  resolve  all  the  scales 
of  motion,  we  have  chosen  N  =  256.  Figures  17  show  the  energy  spectrum 

m- 

Eikj)  =  5:  lu(k,t)l^ 

|k|  =  k 

for  various  times.  It  seems  that  the  dissipation  wavenumber  k,,  is  of  the 
order  of  20.  Hence,  the  previous  estimate  based  on  the  dimension  of  the  at¬ 
tractor  matches  well  the  computational  results.  We  note  also  that,  if  the 
phenomenological  theory  of  turbulence  of  Kraichnan,  predicts  a  decay  of  the 
energy  spectrum  like  k~^,  results  presented  here  seem  to  show  a  faster  decay 
closer  to  k~*.  This  results  are  in  agreement  with  the  ones  obtained  by  Orszag 
in  [23]  and  by  Brachet  et  al.  in  [24],  which  show  that  a  k~^  energy  spectrum 
can  only  be  obtained  when  the  Reynolds  number  is  larger  than  25, 000. 

We  also  remark,  that  if,  at  t  =  0  the  small  scales  corresponding  to  a  wavenum¬ 
ber  larger  than  60  are  set  to  zero,  a  dissipation  range  appears  very  quickly, 
as  we  can  see  at  i  =  5.  The  enstrophy  transfer  from  the  large  scales  to  the 
small  ones  acts  on  the  smallest  scales  after  a  few  iterations.  Figure  18  indi¬ 
cates  that  the  transition  period  is  very  short.  The  small  scales  are  damped 
by  viscous  effect  until  an  equilibrium  between  viscous  and  nonlinear  terms 
appears.  After  that,  we  can  see  that  |  z/y,  I2  oscillates  in  time  and  seems  to 
become  completely  independent  of  its  initial  value.  Figures  28  and  30  repre¬ 
sent  the  isolines  of  the  vorticity,  at  different  times  in  the  interval  [0, 100].  We 
can  see  that  the  very  small  random  structures  of  the  flow  at  the  initial  time 
disappear  quickly.  It  appears  that  fusions  of  these  very  small  structurrs  lead 
to  larger  ones.  So,  after  a  transient  period,  the  flow  is  mainly  constituted  by 
large  structures. 

The  time  step  is  chosen  by  considering  the  accuracy  and  the  stability  of 
the  computation.  For  the  stability.  At  must  satisfy  a  CFL  condition  like 

At  N  \uff  1^00  <  Q  (<  1).  (55) 

From  Figure  25,  we  can  see  that  |  Un  |l«,<  1.25  during  all  the  computation, 
so  (55)  implies  the  following  restriction  on  the  time  step  : 

At  <  1.7  10"^  (a  =  0.5). 
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We  have  set  At  to  10“^.  Here,  the  smallest  scales  (Figure  18)  are  of  the  order 
of  Hence  a  time  step  of  10"^  allows  to  recover  most  of  the  spectrum. 

Indeed,  the  time  differentiation  scheme  used  here  is  a  third  order  method, 
then  the  accuracy  is  of  the  order  of  10~®. 

3.1.2  Comparison  with  a  previous  version  of  the  algorithm. 

In  a  first  time,  we  present  results  obtained  with  a  previous  version  of  the 
multiscale  method.  In  this  version  of  the  algorithm,  was  set  to  (j/  )“* 

and  the  constants  ^i,  $2  of  procedures  (20)  and  (21)  were  a  priori  chosen 
(see  [8]).  This  multiscale  method  is  compared  with  the  pseudo-spectral 
Galerkin  method. 

To  perform  this  analysis,  we  have  retained  two  different  points  of  the  compu¬ 
tational  domain,  namely  x\  =  X2  =  ^  —  l)  and  xj  =  X2  =  ^  —  l) , 

and  we  have  stored  the  value  of  the  horizontal  component  of  the  velocity 
Ui(xi,X2,<)  at  these  points  during  the  time  evolution.  On  Figures  20  and  21, 
we  have  plotted  the  time  history  over  the  interval  [0, 100]  of  those  two  char¬ 
acteristic  values  of  the  flow.  Results  plotted  here  seem  to  be  identical  for 
both  methods,  but  the  differences  between  the  trajectories  obtained  with  the 
different  algorithms  are  too  small  to  appear  on  such  graphic  representation. 
So,  we  have  listed  below  the  exact  values  at  different  intermediate  times  for 
these  trajectories. 

In  Tables  1  and  2,  we  have  listed  the  values  of  the  horizontal  component 
of  the  velocity  at  times  t  =  0,25,50,75  and  100.  These  results  correspond 
to  the  Galerkin  method  for  the  first  column.  In  the  third  column,  we  can 
see  the  difference  between  both  orbits.  It  appears  that  this  quantity  grows 
as  time  evolves  and  becomes  much  larger  than  the  accuracy  which  is  of  the 
order  of  At^  =  10~®  for  this  computation.  On  this  computation,  the  current 
level  Ni{t)  oscillates  between  108  and  200.  Looking  backward  to  the  results 
presented  in  section  2.2  on  the  estimates  of  the  characteristic  length  Tc  (see 
Figures  22),  we  note  that  :  ^ 

j  TNi^  <  10"^  =  At, 

\  <  10-2  ^  lOAt,  for  Ni,  =  108. 

Hence,  levels  Ni^  used  by  the  multiscale  method  are  not  appropriate. 
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Table  1 

Galerkin 

Multiscale 

Difference 

Version  1 

■SH 

mmmmimw 

0.000 

t  =  25 

0.2980770146 

0.2980770543 

o 

II 

0.3592114126 

0.3592114913 

t  =  75 

0.3819454889 

0.3819452151 

<  =  100 

-0.1233715822 

KiiiRiiiriai 

Table  2 

Galerkin 

Multiscale 

Difference 

Version  1 

t  =  0 

0.000 

t  =  25 

0.3879730195 

0.3879729977 

<  =  50 

0.2458823192 

0.2458824078 

lilKlTM 

<  =  75 

-0.3332440190 

-0.3332451088 

<=  100 

0.3196664945 

0.3196686370 

Indeed,  we  recall  that  <  At  means  that,  for  the  level  iVj,  equal  to 
108,  the  scales  smaller  than  can  not  be  fixed  even  on  one  time  iteration. 
Also,TVi,  <  lOAt  means  that  the  coupled  nonlinear  terms  fW*,  ) 

can  not  be  frozen  on  a  time  interval  longer  than  lOAt.  Recalling  now  that 
in  this  version  of  the  multiscale  algorithm,  the  characteristic  length  was  es¬ 
timated  by  therefore  the  length  of  the  frozen  period  oscillated 

between  50  At  and  80  At.  The  constraints  on  the  levels  Nn  and  Ni2  imposed 
by  the  different  estimates  derived  in  section  (2.2)  are  violated  in  this  compu¬ 
tation  and  so,  the  expected  accuracy  can  not  be  recovered  by  the  multiscale 
method. 


Table  3 

Galerkin 

Multiscale 

Difference 

Version  2 

<  =  0 

-0.9650082490  10"^ 

0.0000 

<  =  25 

0.2980770146 

0.2980770146 

<  1.0  10-^° 

<  =  50 

0.3592114126 

0.3592114126 

<  1.0  10-^* 

<  =  75 

0.3819454889 

0.3819454893 

4.0  10-’° 

<  =  100 

-0.1233746011 

-0.1233746004 

7.0  10-** 

50 


Table  4 

Galerkin 

Multiscale 

Difference 

Version  2 

t  =  0 

mmmmmmm 

0.000 

t  =  25 

0.3879730195 

0.3879729977 

t  =  50 

0.2458823192 

0.2458824078 

rfr'ivMM 

f  =  75 

-0.3332440190 

-0.3332451088 

painaiM 

t  =  100 

0.3196664945 

0.319668637U 

3.1..1  Comparison  with  the  improved  version  of  the  algorithm. 

In  Tables  3  and  4,  we  report  results  similar  to  those  in  Tables  1  and  2, 
but  now,  the  second  column  corresponds  to  the  version  of  the  multiscale 
algorithm  presented  in  Section  2.4.  Here,  the  trajectories  of  the  multiscale 
method  remain  close  to  the  trajectories  obtained  with  the  classical  method. 
The  difference  is  less  than  =  10"®  over  the  whole  time  interval  [0, 100]. 
Here,  the  level  iV,i  is  always  larger  than  128,  so  that  the  estimated  limit 
value  of  Tc  is  greater  than  At  (see  Figures  22  and  25).  The  levels  A^,i  and  Ni2 
chosen  by  the  new  algorithm  are  higher  than  those  obtained  by  the  previous 
version.  This  fact  is  due  to  the  restrictions  imposed  by  the  new  criteria  on 
Tc.  For  some  values  of  the  time  t,  the  level  Ni^  decreases  to  a  lower  level 
for  a  short  time  interval  and  then  goes  up  to  its  last  value.  In  such  a  case, 
we  have  remarked  that  the  restriction  imposed  on  Tc  induces  a  restriction  on 
the  level  Ni^.  Indeed,  even  if  a  level  is  acceptable  on  a  few  multigrid  cycles, 
variations  of  the  smallest  scales  or  of  the  transfers  terms  become  too  large, 
so  that  Nil  has  to  change  to  an  upper  level.  Moreover,  we  want  to  mention 
that,  during  the  whole  computation,  neither  the  restriction  due  to  the  small 
scales  evolution  ,  nor  the  restriction  due  to  the  transfer  terms  evolution 
T"Ni^  is  dominant.  Therefore,  it  is  necessary  that  both  criteria  be  retained. 
Finally,  we  want  to  note  that  the  algorithm  can  still  be  improved.  Indeed,  we 
recall  that  the  actual  choice  of  levels  TVj,  and  Ni^  are  determined  according 
two  different  criteria.  One  is  based  on  the  estimate  of  ratios  of  the  kinetic 
energy  (or  enstrophy)  of  the  small  scales  over  the  kinetic  energy  (or  enstro- 
phy)  of  the  large  ones,  and  the  other  one  is  based  on  the  estimates  of  the 
critical  characteristic  time  Tg.  So,  it  follows  that  Tc  may  have  a  relatively 
large  value,  and  then,  the  levels  should  be  adjusted  to  a  lower  value  in  order 
to  have  a  more  reasonable  estimate  of  Tc.  This  can  be  viewed  on  Figure  27, 
where  the  time  evolution  of  Tc  is  plotted.  In  fact,  the  very  strong  oscillations 
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of  Tc  are  not  always  necessary,  and,  by  optimising  the  evaluation  of  Nn  and 
Ni2,  this  time  evolution  can  become  smoother. 

Finally,  to  achieve  the  comparison  between  the  different  algorithms,  we  want 
to  note  the  non-negligible  fact  that  the  multiscale  method  requires  twice 
less  CPU  time  than  the  classical  pseudospectral  method.  On  Figure  24,  the 
quantity  ; 

Tnlg  —  Tg 
Tg 

is  plotted,  where  7/vlg  is  the  CPU  time  required  by  the  nonlinear  Galerkin 
method  (multiscale),  and  Tq  the  CPU  time  required  by  the  classical  Galerkin 
method. 

Finally,  we  want  to  mention  that  these  simulations  required  more  than  75 
hours  of  CPU  on  a  Cray2,  without  counting  all  the  preliminary  tests  needed 
to  the  developments  and  improvements  of  the  algorithm. 
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Figure  19;  Time  evolution  of  |  u^it)  I2. 


Figure  20:  Time  evolution  of  the  horizontal  component  of  the  velocity 
UN(x,y,t)  at  point  x  =  y  =  -  1). 
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Figure  21:  Time  evolution  of  the  horizontal  component  of  the  velocity 
UN{x,y,t)  at  point  x  =  y  =  ^(^  -  1). 


Figure  24:  Time  evolution  of  the  CPU  time  for  the  classical  method  (full 
line)  and  for  the  multilevel  method  (dashed  line). 
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Figure  25:  Time  evolution  of  the  levels  |  ns{t)  |£,oo. 
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Figure  26:  Time  evolution  of  the  levels  Ni^{tj)  and  Ni^{tj) 


Figure  27:  Time  evolution  of  the  characteristic  time  Tc(<j) 


Figure  28:  Vorticity  structures  at  time  t  =  0. 
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Figure  29:  Vorticity  structures  at  time  t  =  30. 
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Figure  30:  Vorticity  structures  at  time  t  =  60. 
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Figme  31:  Vorticity  structures  at  time  t  —  100. 
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3.2  A  direct  numerical  simulation  at  higher  Reynolds 
number 

We  present  here  numerical  results  obtained  in  a  numerical  simulation,  sim¬ 
ilar  to  the  previous  one,  but  with  a  larger  value  of  the  Reynolds  number. 
The  external  force  is  kept  the  same  and  the  viscosity  is  divided  by  4.  The 
initial  condition  is  still  a  random  field  computed  from  a  given  vorticity.  The 
computation  is  performed  over  50,000  iterations,  so  that,  the  flow  is  no  more 
dependant  of  its  unstructured  initial  state.  An  a  posteriori  analysis  of  the 
behavior  of  the  adaptative  multilevel  procedure  confirms  the  previous  as¬ 
sumptions  and  proves  that  the  multiscale  method  is  well  adapted  to  Direct 
Numerical  Simulations.  This  computation  required  about  50  CPU  hours  on 
a  CRAY2.  This  total  computing  time  includes  all  post-treatments  done  by 
the  code,  i.e.  computations  of  different  norms  of  several  quantities  related  to 
the  small  and  large  scales.  The  CPU  time  spent  to  compute  the  velocity  field 
is  32  hours,  which  corresponds  to  7  10“®  second  per  mode  and  per  iteration. 


3.2.1  Description  of  the  initial  condition. 

As  in  (3.1),  the  initial  field  Uo  is  computed  in  the  spectral  space  from  the 
coefficients  of  a  given  vorticity 

i.ij'  =  •£  a„(k)e*",  (56) 

keZ’.jkisAT 


and  the  coefficients  u;o(k)  are  given  by 


wo(k) 


ci9|k|-^e’^»  if  |kl<jb„  =  60 
0  otherwise 


(57) 


where  c\g  is  such  that  |  utg  |^oo  =  2.0.  As  we  can  see  on  Figure  32,  the  slope 
of  the  energy  spectrum  at  the  initial  time  is  equal  to  —3. 

The  viscosity  is  set  to  2.5  10“^,  which  gives  a  Reynolds  number  equal  to 
6, 328  and  the  integral  scale  L  =  2.93.  So,  the  dissipative  wavenumber  kj,  is 
of  the  order  of  28  which  corresponds  to  =  56.  As  we  want  to  p'  ‘brm  a 
Direct  Numerical  Simulation  a  total  number  of  modes  of  512  in  each  dif  iction 
provides  a  grid  fine  enough  to  resolve  the  scales  under  the  dissipative  ones. 
On  Figure  32,  we  can  see  that  the  dissipative  wavenumber  obtained  by  the 
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computation  is  below  60.  Moreover,  on  Figure  32  we  see  that  the  slope  of  the 
energy  si}ectrum  is  smaller  than  —3,  which  is  a  faster  decay  of  the  energy 
than  for  a  fully  develloped  turbulent  6ow. 

As  in  the  previous  computation,  a  time  step  equal  to  10~^  is  small  enough  to 
insure  the  stability  of  the  scheme  and  allows  to  compute  all  the  scales  with 
enough  accuracy. 

3.2.2  Analysis  of  the  computation. 

We  first  want  to  make  some  remarks  on  the  behavior  of  the  multiscale 
method.  Figure  36  shows  the  time  evolution  of  the  two  characteristic  lev¬ 
els  Nil  which  define  the  transition  range.  As  it  was  expected,  the 

variations  of  Ni^  and  follow  the  variations  of  the  ratios 

I  la  .  I  PsxBintiyNi.ZNi)  U 
1  YNi  I2  I  Pni  B{yNi ,  VNt  )  I2 

for  values  of  N\  larger  than  256,  as  we  can  see  on  Figures  33  and  34.  On 
Figures  37  and  38,  the  variations  of  the  characteristic  times  Tsy  and  are 
plotted.  We  remind  that  t^v,  is  used  to  determine  the  lower  level  A/j,  and  the 
length  of  the  period  during  which  the  smallest  scales,  i.e.  for  I  <  fyv.j,  can  be 
frozen.  From  Figure  36,  we  note  that  the  lower  level  of  the  transition  range, 
namely  N,,,  has  to  be  larger  than  256.  In  fact,  we  find  Ni^  of  the  order  of 
320,  for  its  lowest  value.  Figure  35  confirms  that  the  time  derivative  of  the 
velocity  has  a  decaying  spectrum. 

On  Figures  39,  we  have  plotted  the  time  evolution  of  the  different  terms 
appearing  in  the  equation  of  the  small  scale  components  Z/v, .  As  we  have 
observed  in  the  previous  computations,  the  time  derivative  |  za^,  I2  is  of 
the  order  of  the  coupled  nonlinear  terms  |  Qjv,  fiint(yAr,  ,ZAr,)  I2,  while  the 
dissipative  norm  v  ]  Aza^,  I2  is  much  smaller  for  lower  values  of  N\.  For 
the  scales  lying  far  inside  the  dissipation  range,  we  also  note  that  all  the 
quantities  are  of  the  same  order. 

Figure  40  show  the  vorticity  at  different  time  of  the  integration  interval  [0, 50]. 
The  unstructured  initial  field  disappear  completely  after  a  fiew  times  the  eddy 
turnover  time.  After  a  short  transient  period,  the  flow  evolves  by  keeping 
the  same  global  structures.  Figure  40  represents  three-dimensional  views  of 
two-dimensional  maps.  Lighting  technics  have  been  used,  so  that  shadow 
effects  allow  to  see  the  very  fine  structures  of  the  flow. 


62 


Figure  37:  Time  evolution  of  tat,  (c)  = 
iV,  =  480. 
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64,  128,  256  and  M  =  480. 
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Figure  39:  Time  evolution  of  (resp.)  v  \  Az^v,  I2,  |  ijv,  I2,  1  I2 

,  and  I  <3jv,  Biiit(yvv, ,  )  I2  for  Ni  -  64,  128,  256  and  Ni  =  480. 
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3.3  Simulation  of  the  whole  dissipation  range  with 
the  multiscale  method. 

In  this  section,  the  example  is  the  same  as  in  the  Section  3.1.  Our  purpose 
is  to  study  the  effect  of  the  cut-off  value  on  the  computed  solution,  so 
that  several  numerical  experiments  have  been  conducted.  We  have  decreased 
TV,-,  so  that  it  is  of  the  order  of  the  dissipation  wavenumber  k,,  ;  the  whole 
dissipation  range  is  then  modelized  by  the  multiscale  strategy  described  in 
previous  sections.  In  both  simulations,  we  have  noted  that  the  multilevel 
method  allows  to  recover  all  the  large  structures  of  the  flow.  Finally,  we  have 
made  a  similar  test  with  the  classical  method,  i.e.  we  have  decreased  N  in 
order  to  estimate  the  lower  level  required  to  recover  the  large  scales. 


3.3.1  Analysis  of  the  numerical  simulation. 

In  order  to  decrease  the  value  of  TVj,,  the  parameter  e  has  been  set  to  10“*. 
In  this  case,  the  level  TV,,  adjusts  itself,  using  the  procedure  described  in  Sec¬ 
tion  2,  to  the  value  24  which  is  smaller  than  TV,  =  34.  Figure  43  shows  the 
evolution  of  the  two  characteristic  levels  TV,,  and  TV,,,  and  Figure  44  shows 
the  evolution  of  the  characteristic  time  Tc.  The  level  TV,-,  is  quasiconstant  and 
is  equal  to  24  while  the  level  TVj,  is  approximately  equal  to  128  ;  TVj,  is  chosen 
to  be  far  inside  the  dissipation  range.  On  Figures  42,  we  have  represented 
the  energy  spectrum  of  the  computed  solution  at  different  intermediate  times 
of  the  interval  [0,65]  on  which  is  conducted  the  computation.  There  is  no 
energy  pile-up  at  high  wavenumbers  and  the  dissipation  of  the  enstrophy  is 
well  preserved.  By  comparing  the  vorticities  obtained  with  this  simulation 
and  the  Direct  Numerical  Simulation  performed  with  the  classical  Galerkin 
(pseudospectral)  method  with  a  spatial  resolution  equal  to  (256)*,  we  note 
that  all  the  large  scale  vortices  are  well  described  by  the  multilevel  method. 
At  the  times  t  =  60  and  t  =  65,  we  still  have  the  fusion  previously  mentioned. 
We  also  note  that  oscillations  appear  on  the  domain.  They  are  due  to  the 
approximation  made  on  the  small  scales.  They  also  point  out  the  problem 
of  separation  of  scales  related  with  a  Fourier  decomposition.  Indeed,  Fourier 
waves  oscillate  over  the  whole  domain  and  they  can  not  be  directly  associ¬ 
ated  with  a  scale  vorticity.  This  nonlocal  property  implies  the  oscillation 
appearing  on  Figures  47.  Nevertheless,  it  appears  clearly  that  all  the  large 
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Table  7 

uf'^(Xi,X2) 

Difference 

t  =  10 

0.8299845080  10“* 

0.8355520907  10'* 

5.56  10--* 

o 

II 

0.2091267106 

0.2087057712 

4.21  10-“ 

II 

O 

0.3664263833 

0.3652791132 

1.14  10-^ 

t  =  40 

0.2926552919  lO"^ 

0.2587214575  lO"' 

3.39  10-^ 

t  =  50 

0.2458823192 

0.2432537325 

2.63  10-^ 

t  =  60 

0.1573811778 

0.1414418124 

1.59  10-2 

Table  8 

uf^{xi,X2) 

Difference 

t  =  10 

-0.7118241214  10-’ 

-0.7088992727  10"' 

2.92  lO-'* 

t  =  20 

0.1478691156 

0.1470191181 

8.49  lO-'* 

t  =  30 

0.4236670791 

0.4264640060 

2.79  10-^ 

t  =  40 

0.4898684126 

0.4800345247 

9.83  10-^ 

t  =  50 

0.3592114126 

0.3480702322 

1.11  10-^ 

<  =  60 

0.5173406258 

0.5179894219 

6.48  10-'* 

structures  are  captured  with  the  multilevel  method. 

The  total  CPU  time  required  for  this  simulation  is  about  4796  seconds  and 
the  difference  in  (or  energy)  norm  with  the  solution  obtained  with  the 
previous  DNS  is  of  the  order  of  3.6  10“^  at  t  =  50.  The  upper  level  Ni^  of 
the  multilevel  procedure  is  approximatively  equal  to  128  during  the  whole 
computation.  The  small  scale  energy  on  this  level  is  less  than  10~®  and  then 
much  less  than  e.  Hence,  the  level  Ni^  can  be  decre^lsed  to  at  least  64  as  we 
can  see  on  Figure  16  ;  in  that  case,  the  CPU  time  will  be  reduced. 


On  tables  7  and  8,  we  have  compared  the  solutions  obtained  here  with 
the  one  obtained  by  direct  simulation.  It  appears  that  the  difference  between 
the  two  solutions  is  of  the  order  of  10“^  and  hence,  of  the  order  of  e.  The 
oscillations,  appearing  on  Figures  45,  46, and  47  are  probably  due  to  the  fact 
that  the  intermediate  scales,  lying  in  the  inertial  range  and  corresponding 
to  wavenumbers  larger  than  ,  are  not  relaxed  at  each  time  iteration.  A 
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better  understanding  of  the  effects  of  the  multilevel  procedure  on  these  scales 
must  be  done  in  order  to  improve  the  results  presented  here. 


3.3.2  Comparison  with  a  low  resolution  Galerkin  simulation. 

Finally,  in  order  to  show  the  efficiency  of  the  multilevel  method  on  such  sim¬ 
ulations,  we  have  performed  an  additional  test  with  the  classical  method. 
Indeed,  we  have  first  fixed  the  total  number  of  modes  to  32  and  we  have 
integrated  the  system  over  50, 000  time  iterations,  which  corresponds  to  the 
time  interval  [0,50].  The  solution  obtained  is  different  than  the  one  obtained 
when  a  Direct  Simulation  is  done,  i.e.  with  256  modes.  Figure  47  shows 
the  vorticity  structures  obtained  with  this  simulation  at  time  t  =  50.  The 
number  of  modes  is  not  sufficiently  large  so  that  the  large  scales  can  not 
be  computed.  The  vorticity  looks  like  a  non-structured  quantity.  Figure  48 
shows  the  energy  spectrum  at  different  time  of  the  interval  [0,50].  By  fix¬ 
ing  the  resolution  to  (32)^,  it  seems  that  we  do  not  allow  the  appearance  of 
small  scales  and  their  dissipation  mechanism  occuring  in  the  viscous  range, 
the  dynamic  of  the  flow  is  drastically  modified.  We  note  that  this  simulation 
is  stable  in  the  sense  that  no  quantity  grows  artificially.  A  similar  simulation 
has  been  done  with  N  =  4S  instead  of  iV  =  32.  The  picture  of  the  flow  (see 
Figures  47,  49)  is  more  realistic  than  in  the  previous  case.  The  large  scales 
of  the  flow  are  almost  well  captured.  Nevertheless,  the  simulation  performed 
with  the  multilevel  method  with  a  lower  level  A,,  =  24  gives  a  better  quali¬ 
tative  result.  The  CPU  time  required  by  this  simulation  is  2730  seconds  and 
the  difference  with  the  DNS  is  of  the  order  of  3.1  10~^  in  L*  norm.  Hence,  the 
accuracy  recovered  here  is  the  same  as  the  one  obtained  with  the  multilevel 
level  method  when  Ni^  =  24.  The  result  obtained  here  confirms  the  fact  that, 
in  the  previous  computation,  the  level  AT.j  can  be  decreased  to  48  ;  in  such 
a  case,  the  CPU  time  required  by  the  multilevel  method  will  be  at  least  two 
times  less  than  2730  seconds. 

To  conclude,  we  have  proved  with  these  experiences  that  the  multilevel 
method  allows  to  modelize  the  whole  dissipation  range  and  the  end  of  the 
inertial  range  without  disturbing  the  large  scales  of  motion.  Moreover,  we 
have  shown  that  if  the  small  scales  lying  in  the  dissipation  rai  ge  can  be 
modelized,  as  well  as  their  interaction  with  the  large  scales,  they  are  essential 
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to  describe  the  evolution  of  the  flow. 
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Figure  46:  Comparison  of  the  vorticity  structures  at  time  t  =  65. 


Figure  47:  Comparison  of  the  vorticity  structures  obtained  with  the  classical 
Galerkin  method  for  N  =  256  and  N  =  32,  at  time  t  =  50. 
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Figure  49:  Comparison  of  the  vorticity  structures,  between  results  obtained 
with  the  classical  method  for  N  =  256  and  N  =  48,  at  time  t  =  50. 
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Conclusions 


In  this  article,  we  have  proposed  a  multilevel  method  which  treats  differently 
the  large  and  the  small  scales  of  homogeneous  isotropic  flows.  Moreover, 
we  have  derived  mathematical  estimates  for  all  the  parameters  (cut-off  level, 
time  scales)  involved  in  this  dynamical  procedure  leading  to  a  completely 
self-adaptive  procedure.  Firstly,  several  computations  have  been  conducted 
in  the  context  of  DNS,  i.e.  the  whole  spectrum  up  to  the  dissipative  scale  was 
simulated.  In  such  case,  the  multi-level  method  is  able  to  recover  a  velocity 
field  with  the  same  spatial  resolution  than  the  Galerkin  method,  but  with  a 
substantial  speed-up  of  at  least  2  in  CPU  time. 

Secondly,  in  an  approach  similar  to  LES,  we  have  decreased  the  number 
of  modes  retained  for  the  resolved  scales  and  used  the  same  algorithm  to 
modelize  the  interaction  between  the  low  and  high  frequencies.  In  such  case, 
we  have  seen  that  when  the  resolved  scales  are  larger  than  the  dissipative 
ones,  but  of  the  same  order  of  magnitude  (1/12  as  compared  to  1/17),  the 
large  eddies  of  the  flow  are  captured  correctly;  this  is  not  the  case  when  the 
pseudo-spectral  Galerkin  method  is  used  with  only  32^  modes,  i.e.  with  a 
larger  scale  of  the  order  of  1/16. 

Although,  we  only  presented  computational  for  moderately  large  Reynolds 
numbers,  where  the  small  scales  are  strongly  dependent  on  the  energy- 
containing  eddies,  we  think  that  these  results  are  very  promising.  Indeed, 
we  can  reasonably  hope  that  the  multilevel  (nonlinear  Galerkin)  method  can 
be  used  efficiently  for  Large  Eddy  Simulations.  Results  on  this  point  will  be 
presented  elsewhere,  in  the  case  of  two  and  three  dimensional  homogeneous 
isotropic  flows. 

Finally,  we  want  to  observe  that,  by  opposition  with  LES  methods,  no  as¬ 
sumptions  on  the  energy  spectra  or  on  the  velocity  correlations  are  made  here. 
Therefore  the  method  can  be  applied  to  the  simulation  of  non-homogeneous 
flows.  First  attempts  have  been  made  in  this  direction,  for  the  channel  flow 
problem.  Our  efforts  will  be  concentrated  on  this  problem  in  the  near  future. 
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